% set the max. number of refinements:
nrefmax = 6;

% reference solution on finest mesh:
fem_ref = getfem (nrefmax);

% obtain mass matrix for finest mesh:
fema = fem_ref;
fema.equ.c = 0;
fema.equ.a = 1;
fema.xmesh = meshextend (fema);
[Mass,L,M,N] = assemble (fema);
p_ref = fema.mesh.p;
clear fema; % clear to save memory

% interpolate ref. sol. to ref. mesh:
u_ref = postinterp(fem_ref, 'u', p_ref);

% preallocate vectors:
D = repmat (NaN, [nrefmax 1]);
E = repmat (NaN, [nrefmax 1]);

for nref = 0 : nrefmax-1
  fem = getfem (nref);
  D(nref+1) = length(fem.sol.u);

  % interpolate sol. to ref. mesh:
  u_int = postinterp(fem, 'u', p_ref);

  % compute error as a column vector:
  e = u_ref(:) - u_int(:);

  % compute L2-norm of nodal error:
  E(nref+1) = sqrt(e'*Mass*e);
end;

R = E(1:nrefmax-1) ./ E(2:nrefmax);
Q = log2(R);

format compact
format long

D
E
R
Q