output_precision=3;
clear 
x=linspace(0,1,5);
y=linspace(0,2,5);
[nodes,elem,edges]=CreateRectMesh(x,y,1,1,1,1);
[A,b,n2d]=FEMEquationM(nodes,elem,edges,1,0,-1,0,0);
u=FEMSolve(nodes,A,b,n2d,0);
gset zrange [*:*]
gset nokey
ShowSolutionM(nodes,elem,u)
A
b=b'

# evaluate along diagonal
npoints=10;
xv=linspace(0,1,npoints);
yv=linspace(0,2,npoints);
xy=[xv;yv]';

[values, grad] =FEMValueM(xy,nodes,elem,u,0)
