% Finite difference method for solving a boundary value problem of the form v0(x)*d^2y/dx^2+v1(x)dy/dx+v2(x)y=v3(x)
%and the boundary  conditions
%y(x0)=y0 and y(x1)=y1


% solve a boundary value problem of the form v0(x)*d^2y/dx^2+v1(x)dy/dx+v2(x)y=v3(x)
v0 = inline('x^2');
v1 = inline('x');
v2 = inline('x^2');
v3 = inline('0.0');

% Boundary
x0=0.01;
x1=11.7915;
% Boundary conditions
y0=1.0;
y1=0.0;
% Discretize interval [x0,x1]
N=10000;
xx=linspace(x0,x1,N);
dx=(x1-x0)/N;
% initialize linear system resulting from discretization of the BVP
h=zeros(N,N);
y=zeros(N,1);

% fill matrix elements
for i=3:N-2
h(i,i)=-2.0*v0(xx(i))+dx*v1(xx(i))+dx^2*v2(xx(i));
h(i,i-1)=v0(xx(i))-dx*v1(xx(i));
h(i,i+1)=v0(xx(i));
y(i)=dx^2*v3(xx(i));
end

% apply boundary conditions
h(1,1)=1.0;
h(1,2)=0.0;
y(1)=y0;
h(2,2)=-2.0*v0(xx(2))+dx*v1(xx(2))+dx^2*v2(xx(2));
h(2,1)=0.0;
h(2,3)=v0(xx(2));
y(2)=dx^2*v3(xx(2))-(v0(xx(2))-dx*v1(xx(2)))*y0;
h(N-1,N-1)=-2.0*v0(xx(N-1))+dx*v1(xx(N-1))+dx^2*v2(xx(N-1));
h(N-1,N-2)=v0(xx(N-1))-dx*v1(xx(N-1));
h(N-1,N)=0;
y(N-1)=dx^2*v3(xx(N-1))-y1*v0(xx(N-1));
h(N,N)=1.0;
h(N,N-1)=0.0;
y(N)=y1;

%solve linear system
yy=h\y;
%plot solution
plot(xx,yy);