function [T,X] = findiff(p,q,r,a,b,alpha,beta,n)
% [T,X] = findiff(p,q,r,a,b,alpha,beta,n)
% Finite difference solution to the boundary value problem
% x`` = p(t)x`(t)+q(t)x(t)+r(t), x(a) = alpha, x(b) = beta
% This subroutine uses trisys.m to solve a tri-diagonal system.
% p is a function, input.
% q is a function, input.
% r is a function, input.
% a is the left endpoint, input.
% b is the right endpoint, input.
% alpha is the left boundary value, input.
% beta is the right boundary value, input.
% n is the number of steps, input.
% T is the vector of abscissas, output.
% X is the vector of ordinates, output.
T = zeros(1,n+1);
X = zeros(1,n-1);
Va = zeros(1,n-2);
Vb = zeros(1,n-1);
Vc = zeros(1,n-2);
Vd = zeros(1,n-1);
h = (b - a)/n;
for j=1:n-1,
Vt(j) = a + h*j;
end
for j=1:n-1,
Vb(j) = -h^2*feval(r,Vt(j));
end
Vb(1) = Vb(1) + (1 + h/2*feval(p,Vt(1)))*alpha;
Vb(n-1) = Vb(n-1) + (1 - h/2*feval(p,Vt(n-1)))*beta;
for j=1:n-1,
Vd(j) = 2 + h^2*feval(q,Vt(j));
end
for j=1:n-2,
Va(j) = -1 - h/2*feval(p,Vt(j+1));
end
for j=1:n-2,
Vc(j) = -1 + h/2*feval(p,Vt(j));
end
X = trisys(Va,Vd,Vc,Vb);
T = [a,Vt,b];
X = [alpha,X,beta];