function [em,e0]=fdmdirtest16(f,sol,a,b,c,N0,M);
%testing order of convergence of FDM scheme in 1d with Dirichlet bnd i.e.
%-u''+cu=f in (a,b) u(a)= al u(b)=be
%FDM: (xk)_k xk=a+k*h h=(b-a)/N
%(-u(k-1)+2*u(k)-u(k+1))/(h*h)=f(x(k)) k=1..N-1
%%u(0)=al
%u(b)=be
%INPUT
%f,sol - function handles for rhs and solution (we have to know it!)
%a,b - interval ends
%c- constant nonnegative - parameter of PDE
%N0 - size of initial mesh
%M - no of halving step size
%OUTPUT
%e0,em - vectors with discrete L2 and max errors for N=N0,N0*0.5,...,N0*2^(-M+1)
%f=sol=@sin; a=0;b=1;c=0;N0=20;M=6;
%EXAMPLE
%M=6;N0=20;sol=@sin;c=0; a=0;b=1;f=@(x) (1+c)*sin(x);[em,e0]=fdmdirtest16(f,sol,a,b,c,N0,M);
%M=6;N0=20;sol=@(x)(x>0).*(x.*x.*x.*x);c=0; a=-1;b=1;f=@(x)-(x>0).*12.*x.*x;[em,e0]=fdmdirtest16(f,sol,a,b,c,N0,M);
%M=6;N0=20;sol=@(x)(x>0).*(x.*x.*x);c=0; a=-1;b=1;f=@(x)-(x>0).*6.*x;[em,e0]=fdmdirtest16(f,sol,a,b,c,N0,M);
%M=6;N0=20;sol=@(x)(x.*sqrt(x));c=0; a=0;b=1;f=@(x) -0.5*1.5./sqrt(x);[em,e0]=fdmdirtest16(f,sol,a,b,c,N0,M);



 al=sol(a);be=sol(b); 
 em=e0=zeros(M,1);
 N=N0; 
 for k=1:M,
 h=(b-a)/N;
 [u,x]=fdmdir16(f,a,b,al,be,N,c); 
 er=u-sol(x); plot(x,er,";er;"); pause(1);
 em(k)=norm(er,'inf');%max discrete error
 e0(k)=sqrt(h*er'*er); %L2 discrete error
 printf("[%d] %g %g",N,em(k),e0(k));
 if(k>1)
  printf(" \t\t  %f %f\n",em(k-1)/em(k),e0(k-1)/e0(k));
 else printf("\n");
 end%if
 N*=2;
 end%for
end%function


