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