% finite difference solution to a second-order linear ODE 

%** set up model domain
xi=0;			% one end point of domain
xf=1;			% the other end point
N=20;			% number of intervals, N+1 = number of nodes
dx=(xf-xi)/N;	% step size (interval size)
x=[xi:dx:xf]';	% vector holding values of ind. var. at nodes
				% transpose to make a column vector

%* boundary conditions
yi=0;			% y(xi)
yf=0;			% y(xf)

%* dependent variable
y=zeros(N+1,1);
y(1)=yi;			% boundary conditions into y
y(N+1)=yf;


%** build LHS of equation Ay=b,  A matrix

A=zeros(N+1, N+1);

% special cases
A(1,1)= 1;

A(N+1, N+1)= 1;

for i=2:N
   A(i,i-1) = x(i)/dx^2 + 1/dx;
   A(i,i)   = -2*x(i)/dx^2;
   A(i,i+1) = x(i)/dx^2 - 1/dx;
end
   

%** build RHS of Ay=b

b=zeros(N+1,1);		% column vector

% special cases
b(1)= y(1);
b(N+1)= y(N+1);

for i=2:N
  b(i)= -2;
end



%** solve the equation

y=A\b;		


figure(1)
clf
plot(x, y, 'b-'), hold on
title('finite difference solution to x d^2y/dx^2 - 2dy/dx +2 = 0 y(0)=0 y(1)=0')
ylabel('y')
xlabel('x')