page_screen_output=0;
# clear a f zero;

tic
  [nodes,elem,edges]=ReadMesh("./Lshape");
NumberNodes=max(size(nodes))
readingtime=toc
tic
  [A,b,n2d]=FEMEquation(nodes,elem,edges,1,0,'f',0,0);
setuptime=toc
[dof,semiband]=size(A)
tic
  u=FEMSolveSym(nodes,A,b,n2d,0);
solvetime=toc
figure(1);
tic
  ShowSolution(nodes,elem,u)
graphtime=toc
  ShowSolutionMTV(nodes,elem,u,'u.mtv')

gset zrange [*:*]
gset view 60, 300
gset title
gset nokey
replot

print("demo4a.eps","-depsc")

npoints=41;
x=linspace(-1,1,npoints);
y=linspace(1,-1,npoints);
xy=[x;y]';

pause(3)
tic
values=FEMValue(xy,nodes,elem,u,0);
evaluatetime=toc

figure(2);
gset nokey
gset title "A section"
plot (x,values)
print("demo4b.eps","-depsc")

tic
[la,vec]=FEMEig(nodes,elem,edges,'a' ,0,1,4,1e-4);

eigenvaluetime=toc
eigenvalues=la'
figure(3);
gset title "3rd eigenfunction"
ShowSolution(nodes,elem,vec(:,3))
print("demo4d.eps","-depsc")
gset title
# evaluate on a mesh
npoints=21;
xv=linspace(-1,1,npoints);
yv=linspace(-1,1,npoints);

[xx,yy]=meshgrid(xv,yv);
xy=[xx(:),yy(:)];

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

figure(4);
ShowVectorField(xy,-grad)
gset xrange [-1.2:1]
gset yrange [-1:1.2]
replot
print("demo4e.eps","-depsc")
gset xrange [*:*]
gset yrange [*:*]


