function [x, y] = lbvp(p,q,f,n,bcL,bcR)

%  uses centered differences to solve the linear BVP
%       y'' + p(x)y' + q(x)y = f(x),   for xL < x < xR
% where  
%      aL*y(xL) + bL*y'(xL) = cL   and  aR*y(xR) + bR*y'(xR) = cR

% Usage (example at end of file)
%     [x, y] = lbvp(@(x)p(x), @(x)q(x), @(x)f(x),n,bcL,bcR)

% Output: y = solution vector at given x values

% Required arguments
%   p(x): coeff of y' term
%   q(x): coeff of y term
%   f(x): forcing function
%   n: number of x points to be used 
%   bcL = [xL, aL, bL, cL]: values for left BC
%   bcR = [xR, aR, bR, cR]: values for right BC


%%% the code begins
x=linspace(bcL(1),bcR(1),n);
h=x(2)-x(1);
% calculate the coefficients of finite difference equation
for i=2:n-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));
	ff(i)=h*h*f(x(i));
end
if bcL(3)==0
	a(1)=1;  b(1)=0;  c(1)=0;  ff(1)=bcL(4)/bcL(2);
else 
	a(1)=-2+h*h*q(x(1))+2*h*bcL(2)*(1-0.5*h*p(x(1)))/bcL(3);
	b(1)=0;
	c(1)=2;
	ff(1)=h*h*f(x(1))+2*h*bcL(4)*(1-0.5*h*p(x(1)))/bcL(3);
end
if bcR(3)==0
	a(n)=1;  b(n)=0;  c(n)=0;  ff(n)=bcR(4)/bcR(2);
else 
	a(n)=-2+h*h*q(x(n)) - 2*h*bcR(2)*(1+0.5*h*p(x(n)))/bcR(3);
	b(n)=2;
	c(n)=0;
	ff(n)=h*h*f(x(n)) - 2*h*bcR(4)*(1+0.5*h*p(x(n)))/bcR(3);
end

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


% 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


%%% Example: y'' + 100(3x-1)y' + 100xy = 0, for 0 < x < 1, with y(0) = 1, y(1) = 2
%  
% bcL = [0,1,0,1]; bcR = [1,1,0,2];
% [x, y] = lbvp(@(x)100*(3*x-1), @(x)100*x, @(x)0,100,bcL,bcR);
%%% plot solution
% plot(x,y)

% 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  n  (don't be shy about making n 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 has been tested on MATLAB, version R2020b
%  August 8, 2020