%  Solves the steady-state heat equation in a rod with conductivity
%  c(x) = 1 + x^2:
%
%     d/dx( (1+x^2) du/dx ) = f(x),   0 < x < 1
%     u(0) = 1,  u'(1) = 0
%
%  Uses a Galerkin finite element method.

%  Set up grid.

n = input(' Enter number of subintervals: ');
h = 1/n;
x = zeros(n,1);
for i=1:n, x(i)=i*h; end;

%  Form tridiagonal finite element matrix and right-hand side vector.

A = sparse(zeros(n,n));
b = zeros(n,1);

f = inline('2*(3*x^2 - 2*x + 1)');    %  Right-hand side function

%    First equation
x1mh = .5*x(1); x1ph = .5*(x(1)+x(2));
A(1,1) = -(2/h)*(1 + (1/3)*(x(2)^2));
A(1,2) = (1/h)*(1 + (1/3)*(x(2)^2 + x(2)*x(1) + x(1)^2));
b(1) = (h/2)*(f(x1mh)+f(x1ph));      %  Midpoint rule approx to integral
b(1) = b(1) - (1/h)*(1+h^2/3);       %  Contribution of boundary term

%    Equations 2 through n-1
for i=2:n-1,
  ximh = .5*(x(i)+x(i-1));
  xiph = .5*(x(i)+x(i+1));
  A(i,i) = -(2/h)*(1 + (1/3)*(x(i+1)^2 + x(i+1)*x(i-1) + x(i-1)^2));
  A(i,i+1) = (1/h)*(1 + (1/3)*(x(i+1)^2 + x(i+1)*x(i) + x(i)^2));
  A(i,i-1) = A(i-1,i);
  b(i) = (h/2)*(f(ximh)+f(xiph));
end;

%    Last equation
xnmh = .5*(x(n)+x(n-1));
A(n,n) = -(1/h)*(1 + (1/3)*(1 + x(n-1) + x(n-1)^2));
A(n,n-1) = A(n-1,n);
b(n) = (h/2)*(f(xnmh));

%  Solve tridiagonal linear system.

u = A\b;

%  Compare with true solution.  Plot true and computed solution and error.
utrue = (ones(n,1)-x).^2;
err2 = sqrt(h)*norm(u-utrue), % NOTE:  One really should compare at points
errinf = norm(u-utrue,'inf'), % other than the nodes; might get
                              % superconvergence at nodes.
subplot(2,1,1)
plot(x,utrue,'-', x,u,'--')
xlabel('x'); ylabel('u');
title('True solution (solid) and computed solution (dashed)')
subplot(2,1,2)
plot(x,utrue-u,'-')
xlabel('x'); ylabel('Error')