



function [u,x]=fdmdir16(f,a,b,al,be,N,c)
%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 Dirichelt bnd conditions
%c - positive constant parameter of PDE
%N- no od mesh points minus one
%f=@(x)sin(x); c=0; a=0;b=pi; al=be=0; N=20; [u,x]=fdmdir16(f,a,b,al,be,N,c); plot(x,u,x,sin(x),";sol;");
%N=20;sol=@sin;c=3;f=@(x)(1+c)*sin(x); a=0;b=1; al=0;be=sin(b); [u,x]=fdmdir16(f,a,b,al,be,N,c); plot(x,er=u-sin(x),";er;");norm(er,'inf')
%N=20;sol=@(x)(x>0).*(x.*x.*x.*x);c=0; a=-1;b=1;al=sol(a);be=sol(b);f=@(x)-(x>0).*12.*x.*x;[u,x]=fdmdir16(f,a,b,al,be,N,c); plot(x,u,x,sol(x),";sol;");

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)=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=[al;h*h*f(x(2:N));be];

u=A\F;

end

