#!/usr/bin/octave

disp("2D FDM -Laplacian=f on (0,1)^2 ; u=sin(pi*x_1)*sin(pi*x_2) i.e. zero on boundary ");
disp("eh - discrete max norm - testing convergence order in discrete max norm");
disp("el - discrete L^2 norm - testing convergence order in discrete L^2 norm");
clear all
N=7;eh=el=0;
for j=1:5,
 A=Lap2Dsq(N);
 h=1/N;
 M=(N-1)*(N-1);
 f=u=zeros(M,1);

 for k=1:N-1,
  for l=1:N-1,
   u(k+(l-1)*(N-1))=sin(pi*k*h)*sin(pi*l*h);
  endfor
 endfor
 f=h*h*2*pi*pi*u;
 uh=A\f;
 ehp=eh;elp=el;
 eh=norm(u-uh,'inf');
 el=h*norm(u-uh,2);
 printf("N=%4d eh=%6.5e ",N,eh);
 if(j>1)
   printf("e_2h/e_h=%6.5f",ehp/eh); 
 endif
 printf(" el=%6.5e ",el);
 if(j>1)
   printf("el_2h/el_h=%6.5f\n",elp/el); 
 else printf("\n");
 endif
 N*=2;
endfor
