function [em,e0]=fdmdirngtest16(f,sol,a,b,c,N0,M,INCOR);
%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+0.5) so x_N=a+N*h<b - b-x_N=0.5*h
%(-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=1; a=-1;b=1;f=@(x) (1+c)*sin(x);[em,e0]=fdmdirngtest16(f,sol,a,b,c,N0,M);
%INCOR=1;M=6;N0=20;sol=@sin;c=1; a=-1;b=1;f=@(x) (1+c)*sin(x);[em,e0]=fdmdirngtest16(f,sol,a,b,c,N0,M,INCOR);
%M=6;N0=20;sol=@sin;c=1; a=-2;b=1;f=@(x) (1+c)*sin(x);[em,e0]=fdmdirngtest16(f,sol,a,b,c,N0,M);
%INCOR=1;M=6;N0=20;sol=@sin;c=1; a=-2;b=1;f=@(x) (1+c)*sin(x);[em,e0]=fdmdirngtest16(f,sol,a,b,c,N0,M,INCOR);

if(nargin<8)
 INCOR=0;
end
if(INCOR==0)
 printf("RIGHT mesh endpoint not on bnd i.e .x_N<>b - std DIR bnd cond x_N=be\n");
 else
 printf("RIGHT mesh endpoint not on bnd i.e .x_N<>b - Collatz approximation\n");
end

 al=sol(a);be=sol(b); 
 em=e0=zeros(M,1);
 N=N0; 
 for k=1:M,
 h=(b-a)/N;
 [u,x]=fdmdirng16(f,a,b,al,be,N,c,INCOR); 
 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


