

#problem 1
#u''=u u(0)=u(T)=1 by FDM 3 point stencil

bb=[1,2,4,16,20];
for b=bb,
  [u,x]=FDsolver(@(x) 0,a=0,b,al=1,be=1,c=1.0);
  plot(x,u);pause(1);
endfor


#problem 2 test the order of FDM approximation of 
#-u''=sin(x) u(0)=0 u(1)=sin(1) - solution sin(x)
#by doubling the mesh size and computing the ratio e(h)/e(h/2)
#with e(h) the discrete max or L2 norm of the error u_h-r_h u
#plot the error 
N=20;
M=8;
sol=@(x) sin(x); #should be implemented matrix wise i.e. sol(A)=(sol(A(i,j)))
f=@(x) sin(x);
b=1;
erL2=erm=100; 
printf("[N]   erm | erm(N)/erm(2N) | erL2 |  erL2(N)/erL2(2N) \n"); 
for k=1:M,
 [u,x]=FDsolver(f,a=0,b,al=0,be=sol(b),c=0.0,N);
 er=u-sol(x);plot(x,er);pause(0.5);
 perm=erm;
 erm=norm(er,'inf');#max norm
 perL2=erL2;h=(b-a)/N;
 erL2=sqrt(h)*norm(er,2); #discrete L2 norm 
 printf("[%4d] %7.3f ",N,erm);
 if(k>1)  printf("%4.3f ",perm/erm);
 else printf("      ");
 endif
 printf("%7.3f ",erL2);
 if(k>1)  printf("%4.3f \n",perL2/erL2);
 else printf("\n");
 endif
 N*=2;
endfor

#problem 3 
#inexact mesh i.e. x_N=a+N*h<b=x_N+h/2 i.e. with b-x_N<h
#3 poitt stencil on {x_k} approx of -u''
#u(b)=u_N\approx u(x_N)  
#test the order of FDM approximation of 
#-u''=sin(x) u(0)=0 u(1)=sin(1) - solution sin(x)
#by doubling the mesh size and computing the ratio e(h)/e(h/2)
#with e(h) the discrete max the error u_h-r_h u
#plot the error
N=20;
M=8;
sol=@(x) sin(x); #should be implemented matrix wise i.e. sol(A)=(sol(A(i,j)))
f=@(x) sin(x);
b=1;
erL2=erm=100; 
printf("[N]   erm | erm(N)/erm(2N) | erL2 |  erL2(N)/erL2(2N) \n"); 
for k=1:M,
 [u,x]=FDsolver(f,a=0,b,al=0,be=sol(b+(b-a)/(2*N)),c=0.0,N);
 er=u-sol(x);plot(x,er);pause(0.5);
 perm=erm;
 erm=norm(er,'inf');#max norm
 perL2=erL2;h=(b-a)/N;
 erL2=sqrt(h)*norm(er,2); #discrete L2 norm 
 printf("[%4d] %7.3f ",N,erm);
 if(k>1)  printf("%4.3f ",perm/erm);
 else printf("      ");
 endif
 printf("%7.3f ",erL2);
 if(k>1)  printf("%4.3f \n",perL2/erL2);
 else printf("\n");
 endif
 N*=2;
endfor



