function w = linfd ( coeff, a, b, n, alpha, beta )

%LINFD        approximate the solution of the linear boundary value
%             problem
%
%                   u'' = p(x) u' + q(x) u + r(x)
%                   alpha(1) u(a) + alpha(2) u'(a) = alpha(3)
%                   beta(1)  u(b) + beta(2)  u'(b) = beta(3)
%
%             using second-order central finite differences
%
%
%     calling sequences:
%             w = linfd ( coeff, a, b, n, alpha, beta )
%             linfd ( coeff, a, b, n, alpha, beta )
%
%     inputs:
%             coeff   string containing name of m-file defining the 
%                     functions p(x), q(x) and r(x) on the right-hand
%                     side of the differential equation; function should
%                     take a single input and return the values of p, q
%                     and r, in that order; i.e., the m-file header should
%                     be of the form
%                            [p, q, r] = coeff ( x )
%             a       left endpoint of problem domain
%             b       right endpoint of problem domain
%             n       number of uniformly-sized subintervals into which
%                     the domain is to be partitioned
%             alpha   three-component vector of the coefficients which
%                     define the boundary condition at x = a
%             beta    three-component vector of the coefficients which
%                     define the boundary condition at x = b 
%
%     output:
%             w       vector of length n+1 containing the approximate
%                     values of the solution of the boundary value
%                     problem at the locations x = linspace ( a, b, n+1 )
%
%     dependencies:
%             this routine uses the routine TRIDIAGONAL from the "Systems
%             of Equations" library
%

h = ( b - a ) / n;
hsq = h*h;
hd2 = h/2;
x = linspace ( a, b, n+1 );

c   = zeros ( 1, n );
a   = zeros ( 1, n+1 );
b   = zeros ( 1, n );
rhs = zeros ( 1, n+1 );

if ( alpha(2) == 0 )
   a(1) = 1;
   rhs(1) = alpha(3) / alpha(1);
else
    [p q r] = feval ( coeff, x(1) );
	b(1)   = -2;
	rhs(1) = -hsq * r - ( 2 + h * p ) * h * alpha(3) / alpha(2);
	if ( alpha(1) == 0 )
		a(1)   = 2 + hsq * q;
	else
	    a(1)   = 2 + hsq * q - ( 2 + h * p ) * h * alpha(1) / alpha(2);
    end;
end;

for i = 2:n
    [p q r] = feval ( coeff, x(i) );
	c(i-1) = -1 - hd2 * p;
	a(i)   = 2 + hsq * q;
	b(i)   = -1 + hd2 * p;
	rhs(i) = -hsq * r;
end;

if ( beta(2) == 0 )
   a(n+1) = 1;
   rhs(n+1) = beta(3) / beta(1);
else
   [p q r] = feval ( coeff, x(n+1) );
   c(n)   = -2;
   rhs(n+1) = -hsq * r + ( 2 - h * p ) * h * beta(3) / beta(2);
   if ( beta(1) == 0 )
      a(n+1)   = 2 + hsq * q;
   else
      a(n+1)   = 2 + hsq * q + ( 2 - h * p ) * h * beta(1) / beta(2);
   end;
end;

w = tridiagonal ( c, a, b, rhs );
