


function [u,x]=FDMDsolver(f,a,b,al,be,c=0,N=100)
%we want to solve -u''+cu=f in (a,b) with u()=al u(b)=be
%by a FDM method on uniform mesh (x(k)=a+(k-1)*h) k=1,N+1
%u(1)=al bnd cpndition at a
% (1/(h*h)(-u(k-1)+2u(k)-u(k+1)) +c*u(k)=f(x(k)) k=2,..,N  
%u(N+1)=be
%%EXAMPLE
%-u''+u=0 u(0)=u(1)=1
%a=0;b=1;c=1;al=be=1;f=@(x) zeros(size(x));N=100;[u,x]=FDMDsolver(f,a,b,al,be,c)
 h=(b-a)/N;
 x=linspace(a,b,N+1);x=x';
 d=ones(N+1,1);
 A=2*diag(d)-diag(d(1:N),1)-diag(d(1:N),-1);
 A=(1/(h*h))*A;
 A+=c*eye(size(A));
 A(1,2)=A(N+1,N)=0;
 A(1,1)=A(N+1,N+1)=1;

 F=zeros(N+1,1);
 F(1)=al;
 F(2:N)=f(x(2:N));
 F(end)=be;
 u=A\F;
endfunction
