function [em,e0]=fdmmixtest16(f,sol,dsol,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
%(-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,dsol - function handles for rhs and solution and derivative of 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;dsol=@cos;c=0; a=-2;b=1;f=@(x) (1+c)*sin(x);[em,e0]=fdmmixtest16(f,sol,dsol,a,b,c,N0,M);
%INCOR=1;M=6;N0=20;sol=@sin;dsol=@cos;c=0; a=-2;b=1;f=@(x) (1+c)*sin(x);[em,e0]=fdmmixtest16(f,sol,dsol,a,b,c,N0,M,INCOR);

if(nargin<9)
 INCOR=0;
end
if(INCOR==0)
 printf("Neumann bnd at the left end approx. by std diff of order 1\n");
 else
 printf("Neumann bnd at the left end approx. by increased order fdm (using PDE\n");
end
 al=dsol(a);
 be=sol(b);
 em=e0=zeros(M,1);
 N=N0; 
 for k=1:M,
 h=(b-a)/N;
 [u,x]=fdmmix16(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


