
function [u,x,A]=lap1d17(f,a,b,al,be,N,c=1)
%FDM solver for -u''+cu=f in (a,b)
%with u(a)=al u(b)=be (Dirichlet bnd cond)
%INPYT
%f- function handle t ooctave's function computing f(x)
%   (it should work for vectors i.e. f(z)=(f(z(k)) if z a vector
%a,b - ends of interval
%al,be - left/right bnd value
%N - mesh points no minus one
%c - constant coefficient (default =0)
%OUTPUT 
%u - FDM soltion on equidistant mesh
%x- mesh
%EXAMPLE
%f=sol=@sin;c=0;a=0;b=1;al=sol(a);be=sol(b);
%N=20;[u,x]=lap1d17(f,a,b,al,be,N,c);er=u-sol(x);norm(er,'inf')
 x=linspace(a,b,N+1)';
 h=(b-a)/N;
 A=sparse(diag(ones(N-1,1))-diag(ones(N-2,1),1));
 A=A+A'; %interior 3diag matrix with 2 on diag and -1 on sub- super diagonals
 F=(h*h)*f(x(2:N));
 F(1)=F(1)+al;
 F(N-1)=F(N-1)+be;
 B=c*h*h*sparse(eye(N-1))+A;
 u=B\F; %interior values
 u=[al;u;be]; %we area adding bnd values

end