% Usage: 
% y = linfd(p,q,r,a,b,ya,yb,n) or [y x] = linfd(p,q,r,a,b,ya,yb,n)
% Finite difference method for linear boundary value problems of form 
% y"+py'+qy=r
%
% Input:
% p,q,r - Matlab inline functions p(x), q(x), r(x)
% a,b - interval
% ya, yb - boundary conditions
% n - number of panels
%
% Output:
% y - computed solution
% x - mesh points
%
% Examples:
% [y x]=linfd(inline('0'),inline('1'),inline('0'),0,1,1,1,10);
% [y x]=linfd(inline('1/t'),inline('1-1/(4*t^2)'),inline('sqrt(t)*cos(t)'),1,6,1,-0.5,10);
% [y x]=linfd(inline('2'),inline('-3'),inline('t'),-1,1,0,0,10);

function [y x] = linfd(p,q,r,a,b,ya,yb,n)
h = (b - a) / n;
hh = h*h;
h4 = h / 4;
hh2 = hh / 2;
% The first and last rows of the matrix and the right hand side
M = sparse(n+1,n+1);
M(1,1) = 1;
M(n+1,n+1) = 1;
F(1) = ya;
F(n+1) = yb;
x(1) = a;
x(n+1) = b;
% Assemble the matrix and the right hand side
xx = a;
for i = 2 : n
    xx = xx + h;
    Pi = h4 * p(xx);
    M(i,i-1) = -(0.5 - Pi);
    M(i,i) = 1 - hh2 * q(xx);
    M(i,i+1) = -(0.5 + Pi);
    F(i) = hh * r(xx);
    x(i) = xx;
end;
% Solve the resulting linear algebraic equation
y = M \ F';