function phi = basis_eval(Nq,h,x0);

% Evaluate the linear basis functions at the quadrature
% points. 

if Nq==2
  phi = (1/h)*[x0(1); x0(2); 2*h-x0(3); 2*h-x0(4)];
else
  phi = (1/h)*[x0(1); x0(2); x0(3); x0(4); 2*h-x0(5); 2*h-x0(6); ...
               2*h-x0(7); 2*h-x0(8)];
end