
function [x,t]=fdmsolve15(f=@(t) 0,c=@(x) 1.0,a=0,b=1,alpha=1,beta=1,N=100)
#Solves -x''+ cx=f x(a)=alpha x(b)=beta
#Input
#f,c - functions handle
#a,b- the ends of interval
#alpha, beta the bnd values at the ends x(a)=alpha x(b)=beta
#N - the number of subintervals
#Default values f=@(t) 0;c=@(t) 1.0;a=0;b=1;alpha=1;beta=1;N=100;
#Examples
#The simplest call
#x=fdmsolve15; plot(x) #solves -x''+x=0 x(0)=x(1)=1
#[x,t]=fdmsolve15(f=@(t) 0,c=@(x) 1.0,a=0,b=20,alpha=1,beta=1,N=500);plot(t,x)
#solves same ODE but with x(20)=1 - shooting method failed!

h=(b-a)/N;

A1=diag(ones(N-1,1))-diag(ones(N-2,1),1);A1=A1+A1'; 
#A1 tridiagonal with 2 on main diagonal and -1 on sub and super diagnals
A=zeros(N+1,N+1);

A(2:N,2:N)=(1/(h*h))*A1; #puts (1/(h*h))*A1 into A (with zero 1st and last rows/columns)
A(2,1)=A(N,N+1)=-1/(h*h);


#construct rhs and the diagonal c(x_k) component of the matrix 
F=zeros(N+1,1);
F(1)=alpha;F(N+1)=beta;
A(1,1)=A(N+1,N+1)=1.0;
for k=2:N,
 F(k)=f(a+(k-1)*h);
 A(k,k)=A(k,k)+c(a+(k-1)*h);
end

t=linspace(a,b,N+1);
x=A\F;

end


 




