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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function [A,g,xmesh,ucomp,uval] = fe_quadratic(M,a,c,f)
% This function uses a Galerkin finite element method with
% quadratic 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 unknowns, 2 floor(M/2)-1.
%  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 endpoints of the intervals interior to 
%    [0,1] and also the midpoints of the intervals.
%  ucomp is the vector of coefficients of the solution in the
%    finite element basis.
%  uval is the vector of values of the computed u at
%    the points h/2, h, 3h/2, ... 1-h/2, where h = 1/m.
%
% The method uses the weak formulation of the equation
% to compute an approximate solution:
%    a(u,v) = (f,v)  for all piecewise quadratic functions v
% where u(x) = ucomp(1) psi_1(x) + ucomp(2) phi_1(x) + ...
%           + ucomp(2(m-1)-1) psi_{m-1}(x) 
%           + ucomp(2(m-1)) phi_{m-1}(x) + ucomp(2m-1) psi_m(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}].
% The psi basis functions are quadratic hat-functions:
%   psi_j(x_j) = 0, psi_j(x_{j-1}) = 0, psi_j((x_j + x_{j-1})/2) = 1, 
%   and psi_j is zero outside the interval [x_{j-1},x_{j}].
% The number of intervals is m.
% The vector x contains the endpoints of the intervals,
% each of length h.
% fe_quadratic.m     Dianne P. O'Leary   January 2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

m = floor(M/2);
x = linspace(0,1,m+1);
h = x(2);
h2inv = 1/h^2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% The vector xmesh contains the endpoints of the intervals
% interior to [0,1] and also the midpoints of the intervals.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

xmesh = linspace(h/2,1-h/2,2*m-1)';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% We order the basis functions as
%  psi_1, phi_1 , ... psi_{m-1}, phi_{m-1}, psi_m.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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 5 diagonals.
% The (j,k) element is a(v_k,v_j), where v_k is the
%  k-th basis function.
% The j-th element of the right-hand side is (f,v_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 second super diagonal is stored in diag2.
% The data for the right-hand side is stored in rhs.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

diag0 = zeros(2*m-1,1);
diag1 = zeros(2*m-1,1);
diag2 = zeros(2*m-1,1);
rhs   = zeros(2*m-1,1);

for i=1:m
  rhs  (2*i-1) = quad(@fxquad,x(i  ),x(i+1),1.e-6,0, h,x(i)); %(f,psi_i)
  diag0(2*i-1) = quad(@axquad,x(i  ),x(i+1),1.e-6,0, h,x(i)); %a(psi_i,psi_i)

if (i < m)
  rhs  (2*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)); %(f,phi_i)
  diag0(2*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)); %a(phi_i,phi_i)
  diag1(2*i-1) = quad(@axliqu,x(i  ),x(i+1),1.e-6,0, h,x(i)); %a(psi_i,phi_{i})
  diag1(2*i  ) = quad(@axliqu,x(i+1),x(i+2),1.e-6,0,-h,x(i+1)); %a(psi_{i+1},phi_{i})
end
if (i < m-1)
  diag2(2*i  ) = quad(@bxlin,x(i+1),x(i+2),1.e-6,0, h,x(i+1),x(i+2)); %a(phi_i,phi_{i+1})
end
end

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

diag1u = [0;  diag1(1:2*m-2)];
diag2u = [0;0;diag2(1:2*m-3)];
A = spdiags([diag2,diag1,diag0,diag1u,diag2u],-2:2,2*m-1,2*m-1);
g = rhs;

ucomp = A \ g;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Compute the function value at the m-1 interior mesh points and
% the midpoints of each of the intervals.
% The value at the midpoint of an interval is the coefficient of
% the quadratic basis function plus the values contributed by the
% one or two linear basis functions that are nonzero there.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

uval= ucomp;
uval(1) = uval(1) + .5*uval(2);

for i=3:2:2*m-3,
  uval(i) = uval(i) + .5*(uval(i-1)+uval(i+1));
end

uval(2*m-1) = uval(2*m-1) + .5*uval(2*m-2);