



function [u,x,A,F]=fdmmix16(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
%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 Neumann/Dirichlet bnd conditions
%c - positive constant parameter of PDE
%N- no od mesh points minus one
%INCOR- 0 standard forward diff approx of u'(a) or increased order by using PDE
%     INCOR=0 -> (u(1)-u(0))/h=al; INCOR=1 -> (u(1)-u(0))/h - 0.5*h*c*u(0)=al - 0.5*h*f(a)
%                                           using (u(a+h)-u(a))/h=al + 0.5*h*u''(a)+O(h^2)=
%                                                                =al+0.5*h*(c*u(a)-f(a))   
%EXAMPLE
%s=@(x)sin(x);ds=@(x)cos(x);c=0;f=@(x)(1+c)*s(x);a=0;b=pi;al=ds(a);be=s(b);N=20;[u,x]=fdmmix16(f,a,b,al,be,N,c); plot(x,u-sin(x),";er;");
%N=20;s=@(x)sin(x);ds=@(x)cos(x);c=1;f=@(x)(1+c)*s(x); a=-1;b=2; al=ds(a);be=s(b); [u,x]=fdmmix16(f,a,b,al,be,N,c); plot(x,er=u-s(x),";er;");norm(er,'inf')%repeat for N=40
%INCOR=1.0;N=20;s=@(x)sin(x);ds=@(x)cos(x);c=1;f=@(x)(1+c)*s(x); a=-1;b=2; al=ds(a);be=s(b); [u,x]=fdmmix16(f,a,b,al,be,N,c,INCOR); plot(x,er=u-s(x),";er;");norm(er,'inf')%repeat for N=40


if(nargin<8)
 INCOR=0;
end;

h=(b-a)/N;
%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)) ;

A=sparse(zeros(N+1,N+1));
A(1,1)=-1.0  -  0.5*INCOR*c*h*h; 
A(1,2)=1.0;
A(N+1,N+1)=1;
A(2,1)=A(N,N+1)=-1;
A(2:N,2:N)=iA + h*h*c*speye(N-1,N-1);

%mesh
x=linspace(a,b,N+1); x=x';

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

u=A\F;

end

