By Mark S. Gockenbach (siam, 2010)

K=Stiffness1(T,a);

F=Load1(T,f);

U=K\F;

V=NodalValues1(T,u);

ShowPWLinFcn1(T,V-U)

K=Stiffness1(T,a);

F=Load1(T,f);

U=K\F;

V=NodalValues1(T,u);

ShowPWLinFcn1(T,V-U)

Using the code

1. 
   Write a command called Mass1 that computes the mass matrix. The calling sequence should be simply M=Mass1(T), where T is the triangulation.
   
2. 
   Choose some other geometric shapes and/or combinations of boundary conditions, and write some mesh generation routines analogous to RectangleMeshD1.
   
3. 
   Extend the code to handle inhomogeneous Dirichlet conditions. Recall that such boundary conditions change the load vector, so the routine Load1 must be modified.
   
4. 
   Extend the code to handle inhomogeneous Neumann conditions. Like inhomogeneous Dirichlet conditions, the load vector is affected
   
5. 
   (Hard) Write a routine to refine a given mesh, according to the standard method suggested in Exercise 10.1.4 of the text.
   

As I mentioned above, the mesh data structure described in Mesh1.m includes the information needed to solve exercises 3, 4, and 5.