function lbvp

%  uses centered differences to solve the linear BVP
%       y'' + p(x)y' + q(x)y = f(x)   for xL < x < xr
% where  
%      a0*y(xl) + b0*y'(xl) = c0   and  a1*y(xR) + b1*y'(xR) = c1

% it is required that a0 and b0 are not both zero (and same for a1 and b1)

% an explanation of how to use this code is given at the end of this file


% set boundary conditions
	xL=0; a0=1; b0=0; c0=1;
	xR=1; a1=1; b1=0; c1=2;

% specify number of grid points
nx=100;


x=linspace(xL,xR,nx);
h=x(2)-x(1);
% calculate the coefficients of finite difference equation
for i=2:nx-1
	a(i)=-2+h*h*q(x(i));
	b(i)=1-0.5*h*p(x(i));
	c(i)=1+0.5*h*p(x(i));
	f(i)=h*h*rhs(x(i));
end;
if b0==0
	a(1)=1;  b(1)=0;  c(1)=0;  f(1)=c0/a0;
else 
	a(1)=-2+h*h*q(x(1))+2*h*a0*(1-0.5*h*p(x(1)))/b0;
	b(1)=0;
	c(1)=2;
	f(1)=h*h*rhs(x(1))+2*h*c0*(1-0.5*h*p(x(1)))/b0;
end;
if b1==0
	a(nx)=1;  b(nx)=0;  c(nx)=0;  f(nx)=c1/a1;
else 
	a(nx)=-2+h*h*q(x(nx)) - 2*h*a1*(1+0.5*h*p(x(nx)))/b1;
	b(nx)=2;
	c(nx)=0;
	f(nx)=h*h*rhs(x(nx)) - 2*h*c1*(1+0.5*h*p(x(nx)))/b1;
end;


% solve the tri-diagonal matrix problem
y=tri(a,b,c,f);

% plot solution
plot(x,y)
hold on
box on
grid on
%axis([0 1.02 -2 2])

xlabel('x-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')
set(gca,'FontSize',14);


% define functions in BVP

function g=p(x)
g=100*(3*x-1);

function g=q(x)
g=100*x;

function g=rhs(x)
g=0;



% tridiagonal solver
function y = tri( a, b, c, f )
N = length(f);
v = zeros(1,N);   
y = v;
w = a(1);
y(1) = f(1)/w;
for i=2:N
    v(i-1) = c(i-1)/w;
    w = a(i) - b(i)*v(i-1);
    y(i) = ( f(i) - b(i)*y(i-1) )/w;
end
for j=N-1:-1:1
   y(j) = y(j) - v(j)*y(j+1);
end




%  This code is run by simply entering the command  lbvp  in MATLAB (making 
%  sure that this file is in the path for MATLAB).  The user needs to first 
%  enter certain information about the problem into this file before running it.

%  Needed Input From User

%  Boundary conditions (BCs):  xL, a0, b0, c0 (line 14)  and xR, a1, b1, c1 (line 15)
%  Number of grid points:  nx (line 18)
%    
%  p(x):  this is the coefficient of y' (line 66)
%  q(x):  this is the coefficient of y (line 69)
%  rhs(x):  this is the right hand side of the equation (line 72)

%  Example:   y'' + 100(3x-1)y' + 100xy = 0  for 0 < x < 1  with y(0) = 1 and y(1) = 2
%    xL=0; a0=1; b0=0; c0=1;
%    xR=1; a1=1; b1=0; c1=2;
%    nx = 100
%    p = 100*(3x-1)
%    q = 100*x 
%    rhs = 0

% What to do if it doesn't work (aside from checking on input errors)
%  1.  If the problem has a boundary or interior layer you should consider 
%      increasing  nx  (don't be shy about making nx very large).  If the layer
%      is sharp then you might want to try  sbvp.m  
%  2.  If q(x) is negative in the interval, then the possibility arises that 
%      the solution is either non-unique or non-existent

%  Background
%  This is a simple code that solves linear 2nd order BVPs. It does not 
%  use adaptive methods, but instead uses 2nd order centered differences on 
%  a uniform grid.  The specific algorithm is described in the book Intro to Numerical  
%  Methods for Differential Equations (Holmes, 2007) on pages 46-49.  It is very  
%  robust, and was used to solve all of the linear BVPs in Intro to Perturbation 
%  Methods (Holmes, 1995).  It is also relatively fast, and as an illustration 
%  is solves the above example in 0.004 sec on a laptop.  If you replace the 
%  100's in the DE with 10000's, and take nx = 10000, it takes only 0.48 sec. 

%  It has been tested on MATLAB, version R2009b
%  April 8, 2010