function [A,g,xmesh,ucomp] = finitediff2(M,a,c,f)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [A,g,xmesh,ucomp] = finitediff2(M,a,c,f)
% This function uses a 2nd order finite difference method
% to form an approximate solution to the differential equation
%   -(a(x) u'(x))' + c(x) u(x) = f(x)  for x in (0,1)
% with u(0)=u(1)=0.
%
% The function a should satisfy a(x) \ge a_0 > 0 on (0,1),
% and the function c should satisfy c(x) \ge 0 on (0,1).
%
% Inputs:
%  M  defines the number of mesh points, including the
%     endpoints of the interval.
%  a, c, and f are user-defined Matlab functions.
%   Their input is a row or column vector of x coordinates
%   in the interval [0,1], and their output is a vector of 
%   values of a(x), c(x), or f(x).  
%   The function a also has a second output argument, a vector of 
%    values of a'(x).
%
% Outputs:
%  A is the matrix from which the coefficients of the solution are
%    computed.
%    The matrix entries are computed using central difference 
%    approximations for u'' and u':
%          u'(x) -> (u(x+h)-u(x-h)) / (2h)
%       - u''(x) -> (-u(x+h)+ 2 u(x) -u(x-h)) / (h^2)
%  g is the right-hand side of the matrix problem.
%  xmesh contains the equally-spaced points at which the solution 
%    is computed.
%  ucomp is the vector of values of the
%    computed solution at the interior mesh points.
%
% finitediff2.m     Dianne P. O'Leary   January 2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize some parameters.
%   h is the spacing between meshpoints.
%   x is the vector of meshpoints.
%   xmesh is the vector of interior meshpoints.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = linspace(0,1,M);
h = x(2);

h2inv = 1/h^2;
n = M - 2;

xmesh = x(2:n+1)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Form the diagonals of A 
%   and form the right-hand side rhs so that 
%            A ucomp = rhs.
% A is a tridiagonal matrix with values 2a/h^2 + c on its main diagonal,
% -a/h^2 + a'/(2h) on its subdiagonal,
% and -a/h^2 - a'/(2h) on its superdiagonal.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[a0,ap] = feval(a,xmesh);
adiff = a0*h2inv;
bdiff = ap/(2*h);
cdiff = feval(c,xmesh);
rhs   = feval(f,xmesh);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Pad the vectors to make Matlab's spdiags pick up the proper
% values for each row.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

adiffu = [0;adiff(1:n-1)];
bdiffu = [0;bdiff(1:n-1)];
bdiffl = [bdiff(2:n);0];
adiffl = [adiff(2:n);0];
A = spdiags([-adiffl+bdiffl,2*adiff + cdiff,-bdiffu-adiffu],-1:1,n,n);
g = rhs;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Solve the linear system.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

ucomp = A \ g;