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];
