  % ... set up boundary condition matrix:
  %     homogeneous Dirichlet conditions: h u = r with h = 1 and r = 0.
  v = [ 1 ; ... % = scalar equation
        1 ; ... % = number of Dirichlet conditions
        1 ; ... % = length of q
        1 ; ... % = length of g
        1 ; ... % = length of h
        1 ; ... % = length of r
       '0'; ... % = value of q
       '0'; ... % = value of g
       '1'; ... % = value of h
       '0'];    % = value of r
  b = [v v v v]; % Dirichlet conditions on all four boundary segments

  % ... set up point matrix p, edge matrix e, and triangle matrix t:
  p = [1 2 0 1 2 3 3 3 2 1 0 0; ... % x-coordinates of vertices
       1 1 0 0 0 0 1 2 2 2 2 1];    % y-coordinates of vertices
  % for each boundary segment, we have one column with ...
  % (Note: subdomain = 0 indicates outside boundary.)
  e = [ 3   4   5   6   7   8   9  10  11  12 ; ... % starting vertex
        4   5   6   7   8   9  10  11  12   3 ; ... % ending vertex
        0  1/3 2/3  0  1/2  0  1/3 2/3  0  1/2; ... % starting parameter
       1/3 2/3  1  1/2  1  1/3 2/3  1  1/2  1 ; ... % ending parameter
        1   1   1   2   2   3   3   3   4   4 ; ... % boundary segment number
        1   1   1   1   1   1   1   1   1   1 ; ... % left subdomain
        0   0   0   0   0   0   0   0   0   0 ];    % right subdomain
  % (Note: vertices must be in counter-clockwise ordering.)
  t = [  1  1  1  2  2  5 10  1  1  1  2  2; ... % local vertex 1
        12  3  4  4  5  6 11 10  9  2  8  7; ... % local vertex 2
         3  4  2  5  7  7 12 12 10  9  9  8; ... % local vertex 3
         1  1  1  1  1  1  1  1  1  1  1  1];    % subdomain number

  % ... test mesh by printing out element and vertex numbers:
  set (0, 'DefaultAxesFontSize',16);
  plotmeshtest
  set (gca, 'xtick',[0:1/3:3], 'ytick',[0:1/3:2]);
  title (['Grid for small test case']);
  set (0, 'DefaultAxesFontSize','factory');

  % ... set problem coefficients:
  %     - div (c grad(u)) + a u = f
  c = 1;
  a = 0;
  f = 1;

  % ... assemble stiffness matrix and right-hand side vector:
  %     Note: This assembly has eliminated the Dirichlet boundary conditions
  %     from the linear system of equations.
  [K, F, B, ud] = assempde (b, p,e,t, c,a,f);

  % ... solve system and re-assemble solution vector on all vertices:
  u1 = K \ F; % solve linear system for solution on interior vertices only
  u = B * u1 + ud; % re-assemble solution on all vertices

  % ... surface plot with lines for triangle edges:
  %     Note: For finer meshes, it becomes impossible to show all
  %     edges, but for smaller meshes it is quite instructive.
  plotsolution
  title (['Solution to small test case']);

