



function [u,x]=fdmdirng16(f,a,b,al,be,N,c,INCOR)
%fdm solver for -u''+cu=f in (a,b) (c>=0 nonnegative constant)
%u(a)=al u(b)=be
%using nonconforming grid i.e .h=(b-a)/(N+0.5) so x_N<b i.e. b-x_N=0.5*h
%INPUT
%f - function handle to rhs f (it must work for vectors i.e .f(x)=[f(x(1),...,f(x(N-1))]
%a,b- interval ends
%al,be - rhs of Dirichelt bnd conditions
%c - positive constant parameter of PDE
%N- no od mesh points minus one 
%INCOR - =1  increased order at the right end by Collatza approximation
%        =0 standard Dirichlet bnd i.e. x_N=be (x_N \approx u(x_N) 
%EXAMPLE
%N=20;s=@(x)sin(x); c=0;f=@(x)(1+c)*s(x); a=0;b=pi; al=s(a);be=s(b);[u,x]=fdmdirng16(f,a,b,al,be,N,c); plot(x,er=u-s(x),";sol;");norm(er,'inf')
%INCOR=1;N=20;s=@(x)sin(x); c=0;f=@(x)(1+c)*s(x); a=0;b=pi; al=s(a);be=s(b);[u,x]=fdmdirng16(f,a,b,al,be,N,c,INCOR); plot(x,er=u-s(x),";sol;");norm(er,'inf')
%N=20;sol=@sin;c=3;f=@(x)(1+c)*sin(x); a=-2;b=0.5; al=sol(a);be=sol(b); [u,x]=fdmdirng16(f,a,b,al,be,N,c); plot(x,er=u-sin(x),";er;");norm(er,'inf')
%INCOR=1;N=20;sol=@sin;c=3;f=@(x)(1+c)*sin(x); a=-2;b=0.5; al=sol(a);be=sol(b); [u,x]=fdmdirng16(f,a,b,al,be,N,c,INCOR); plot(x,er=u-sin(x),";er;");norm(er,'inf')

if(nargin<8)
 INCOR=0;
end

h=(b-a)/(N+0.5);
%3 diagonal matrix N-1 X N-1 with 2 on main diag and -1 on sub super diags
iA=2*sparse(diag(ones(N-1,1))) - sparse(diag(ones(N-2,1),1)) - sparse(diag(ones(N-2,1),-1)) ;

%mesh
x=a+h*(0:N); x=x';

A=sparse(zeros(N+1,N+1));
A(1,1)=1.0;
if(INCOR==0)
 A(N+1,N+1)=1;
else
 %l(t)=u(N+1)*(t-x(N))/h + u(N)*(x(N+1)-t)/h and the last equation l(b)=be
 %u(N+1)*(b-x(N))/h + u(N)*(x(N+1)-b)/h=u(N+1)*1.5 - u(N)*0.5 = be
 A(N+1,N+1)=(b - x(N))/h;
 A(N+1,N)  =(x(N+1)-b)/h;
 %testing%B=full(A);B(N+1,N-1:N+1)
end

A(2,1)=A(N,N+1)=-1;
A(2:N,2:N)=iA + h*h*c*speye(N-1,N-1);



%h scaling goes also into rhs .e. F
F=zeros(N+1,1);
F=[al;h*h*f(x(2:N));be];

u=A\F;

end

