function z = axliqu(xpts,h,xend)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% function z = axliqu(xpts,h,xend)
% This function computes the values of 
%    a phi'_i psi'_i + c phi_i  psi_i          if h>0
%    a phi'_{i-1} psi'_i + c phi_{i-1} psi_i   if h<0
% at the points specified by xpts,
% where phi is a linear basis element and psi is a quadratic 
% basis element for the finite element method.
% It is used by quad in computing matrix entries for the finite element
% formulations.
%
% Input:
%  xpts is a row or column vector of points
%  abs(h) is the distance between mesh points 
%  xend is an endpoint of the interval over which the basis function
%    is nonzero.
%
% The definition of phi is
%             phi(x) = (x - xend)/h,
% so if h>0, then xend is the right endpoint of the interval over
% which the basis function is nonzero, and we are looking at the
% interval over which phi is increasing.
% If h < 0, then xend is the left endpoint, and we are looking at the
% interval over which phi is decreasing.
% The definition of psi is
%             psi(x) = -4 (x - xend)(x - (xend+h))/h^2.
%
% Functions required: a, c
% Outputs:
%  z =  a(xpts) phi'_i(xpts) psi'_i(xpts) 
%        + c(xpts) phi_i(xpts)  psi_i(xpts)          if h>0
%  z =  a(xpts) phi'_{i-1}(xpts) psi'_i 
%        + c(xpts) phi_{i-1}(xpts) psi_i(xpts)       if h<0
%
% axliqu.m      Dianne P. O'Leary   January 2005
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

habs = abs(h);
psi      = -4*(xpts-xend).*(xpts-(xend+habs))/h^2;
psiprime = -4*(2*xpts-2*xend-habs)/h^2;

if (h < 0)
  phi =  (xpts-xend-habs)/h;
else
  phi =  (xpts-xend)/h;
end

phiprime = 1/h;

z = a(xpts).*psiprime*phiprime + c(xpts).*psi.*phi;