a=0;b=1;c=0;d=1;
f=inline('8*pi^2*sin(2*pi*x).*cos(2*pi*y)','x','y');
g=inline('sin(2*pi*x).*cos(2*pi*y)','x','y');
uex=g; nx=4; ny=4;
for n=1:5
  [u,x,y,error(n)]=poissonfd(a,c,b,d,nx,ny,f,g,uex);
  nx = 2*nx; ny = 2*ny;
end
format short e; error
assert(error,[1.3565e-01   4.3393e-02   1.2308e-02   3.2775e-03   8.4557e-04],5e-6);
p=log(abs(error(1:end-1)./error(2:end)))/log(2)
assert(p,[1.6443e+00   1.8179e+00   1.9089e+00   1.9546e+00],1e-4);