function [u,utrue,errL2,errH1] = fem1d()
%
%  program to solve the equation
%
%     -d/dx(p(x)du/dx) + q(x)u = f
%
%  with boundary conditions
%
%     u(0) = u(pi) = 0
%
%  using the finite element method with piecewise linear
%  approximation functions.   Uses functions p, q, and f
%  to evaluate coefficients and right-hand side.
%
%  User is asked for the number of subintervals n.
%  Routine plots the true and computed solutions and their
%  difference.  Returns the computed solution u, the true
%  solution utrue, the L2-norm of the error errL2, and the
%  H1-seminorm of the error errH1.
%
%  Stiffness matrix is stored as a dense matrix and solved
%  using MATLAB operator "\".
%

n = input(' Enter number of subintervals: ');
h = pi/n;

%  Loop over subintervals.  Assemble global stiffness matrix A
%  and right-hand side vector b.

A = zeros(n-1,n-1);  b=zeros(n-1,1);
for i=1:n,
  xim1 = (i-1)*h;
  xi = i*h;
  pt1 = (xi+xim1)/2 - h/(2*sqrt(3));  % pt1 and pt2 are the points used
  pt2 = (xi+xim1)/2 + h/(2*sqrt(3));  % in 2-point Gauss quadrature for
                                      % evaluating integrals
%   Stiffness matrix
  pint = .5*(p(pt1)+p(pt2))*h;
  qiint = .5*(q(pt1)*(pt1-xim1)^2 + q(pt2)*(pt2-xim1)^2)*h;
  qim1int = .5*(q(pt1)*(xi-pt1)^2 + q(pt2)*(xi-pt2)^2)*h;
  if i < n, A(i,i) = A(i,i) + pint/h^2 + qiint/h^2; end;
  if i > 1,
    A(i-1,i-1) = A(i-1,i-1) + pint/h^2 + qim1int/h^2;
    qiim1 = .5*h*(q(pt1)*(xi-pt1)*(pt1-xim1) + q(pt2)*(xi-pt2)*(pt2-xim1));
    if i < n,
      A(i-1,i) = -pint/h^2 + qiim1/h^2;
      A(i,i-1) = A(i-1,i);
    end;
  end;

%   Right-hand side
  fiint = .5*(f(pt1)*(pt1-xim1) + f(pt2)*(pt2-xim1))*h;
  fim1int = .5*(f(pt1)*(xi-pt1) + f(pt2)*(xi-pt2))*h;
  if i < n, b(i) = b(i) + fiint/h; end;
  if i > 1, b(i-1) = b(i-1) + fim1int/h; end;
end;

%  Solve linear system.

u = A\b;

%  Plot true and approximate solutions and error.

for i=1:n-1,
  xi = i*h;
  if xi <= pi/3, utrue(i,1) = sin(3*xi); end;    
  if xi > pi/3, utrue(i,1) = .01*sin(3*xi); end;
end;
plot([h:h:pi-h]',u,'--', [h:h:pi-h]',utrue,'-')
xlabel('x')
ylabel('u(x)')
title('True solution (solid) and approximate solution (dashed)')
pause(5),
plot([h:h:pi-h]',utrue-u,'-')
xlabel('x')
ylabel('error')
title('Difference between true and computed solution')


%  Compute (approximate) L2-norm and H1-seminorm of error.

errL2 = 0;  errH1 = 0;
for i=1:n,
  xim1 = (i-1)*h;
  xi = i*h;

  pt1 = (xi+xim1)/2 - h/(2*sqrt(3));  % pt1 and pt2 are the points used
  pt2 = (xi+xim1)/2 + h/(2*sqrt(3));  % in 2-point Gauss quadrature for
                                      % computing the error.
  if i > 1 & i < n,
    u_pt1 = u(i-1)*(xi-pt1)/h + u(i)*(pt1-xim1)/h;  % Compute function
    u_pt2 = u(i-1)*(xi-pt2)/h + u(i)*(pt2-xim1)/h;  % values and derivs
    uprime = (u(i)-u(i-1))/h;                       % for piecewise linear
  elseif i==1,                                      % approximate solution
    u_pt1 = u(i)*(pt1-xim1)/h;                      % at pt1 and pt2.
    u_pt2 = u(i)*(pt2-xim1)/h;
    uprime = u(i)/h;
  elseif i==n,
    u_pt1 = u(i-1)*(xi-pt1)/h;
    u_pt2 = u(i-1)*(xi-pt2)/h;
    uprime = -u(i-1)/h;
  end;
  if pt1 <= pi/3,                                   % Compute function
    utrue_pt1 = sin(3*pt1);                         % values and derivs
    utrueprime_pt1 = 3*cos(3*pt1);                  % for true solution
  else                                              % at pt1 and pt2.
    utrue_pt1 = .01*sin(3*pt1);
    utrueprime_pt1 = .03*cos(3*pt1);
  end;
  if pt2 <= pi/3,                                   
    utrue_pt2 = sin(3*pt2);                        
    utrueprime_pt2 = 3*cos(3*pt2);                  
  else                                             
    utrue_pt2 = .01*sin(3*pt2);
    utrueprime_pt2 = .03*cos(3*pt2);
  end;

  errL2 = errL2 + h*.5*( (u_pt1 - utrue_pt1)^2 + (u_pt2 - utrue_pt2)^2 );
  errH1 = errH1 + h*.5*( (uprime - utrueprime_pt1)^2 + ...
                         (uprime - utrueprime_pt2)^2 );
end;

errL2 = sqrt(errL2),  errH1 = sqrt(errH1),

%
%  Functions defining coefficients and right-hand side
%
function pval = p(x)
if x <= pi/3,
  pval = 1;
else
  pval = 100;
end;
return;

function qval = q(x)
qval = 1;
return;

function fval = f(x)
sin3x = sin(3*x);
fval = 9*sin3x;
if x <= pi/3,
  fval = fval + sin3x;
else
  fval = fval + .01*sin3x;
end;