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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [A,g,xmesh,ucomp] = fe_linear(M,a,c,f)
% This function uses a Galerkin finite element method with
% linear elements 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.
%  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 coefficients of the solution in the
%    finite element basis, equal to the value of the
%    computed solution at the interior mesh points.
%
% The method uses the weak formulation of the equation
% to compute an approximate solution:
%    a(u,v) = (f,v)  for all piecewise linear functions v
% where u(x) = ucomp(1) phi_1(x) + ... + ucomp(M-2) phi_{M-2}(x)
% Here,
%  a(u,v) = int(a(x) u'(x) v'(x) + c(x) u(x) v(x)) dx,
%   (f,v) = int(f(x) v(x)) dx,
% and the limits of integration are [0,1].
%
% The phi basis functions are linear hat-functions:
%   phi_j (x_j) = 1, phi_j (x_k) = 0 for j not equal to k,
%   and phi_j is zero outside the interval [x_{j-1},x_{j+1}].
% fe_linear.m     Dianne P. O'Leary   January 2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Initialize some parameters.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

x = linspace(0,1,M);
h = x(2);
n = M - 2; % n is the number of unknowns.
xmesh = x(2:n+1)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the matrix and the right-hand side in order
%  to solve for the coefficients ucomp of the approximate solution.
% The matrix is symmetric, made up of 3 diagonals.
% The (j,k) element is a(phi_k,phi_j).
% The j-th element of the right-hand side is (f,phi_j).
%
% The data for the main diagonal is stored in diag0.
% The data for the first superdiagonal is stored in diag1.
% The data for the right-hand side is stored in rhs.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

diag0 = zeros(n,1);
diag1 = zeros(n,1);
rhs   = zeros(n,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Form the matrix and solve the linear system.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for i=1:n,
diag0(i) = quad(@axlin,x(i  ),x(i+1),1.e-6,0, h,x(i)) ...
      + quad(@axlin,x(i+1),x(i+2),1.e-6,0,-h,x(i+2));
diag1(i) = quad(@bxlin,x(i+1),x(i+2),1.e-6,0, h,x(i+1),x(i+2));
rhs  (i) = quad(@fxlin,x(i  ),x(i+1),1.e-6,0, h,x(i  )) ...
      + quad(@fxlin,x(i+1),x(i+2),1.e-6,0,-h,x(i+2));
end

diag1u= [0;diag1(1:n-1)];
A = spdiags([diag1,diag0,diag1u],-1:1,n,n);
g = rhs;

ucomp = A \ g;