%  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) = 0,  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,1);
for i=1:n-1, x(i)=i*h; end;

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

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

f = inline('-2*(3*t^2 - t + 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

%    Equations 2 through n-2
for i=2:n-2,
  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
xnm1mh = .5*(x(n-1)+x(n-2));
xnm1ph = .5*(x(n-1)+1);
A(n-1,n-1) = (2/h)*(1 + (1/3)*(1 + x(n-2) + x(n-2)^2));
A(n-1,n-2) = A(n-2,n-1);
b(n-1) = -(h/2)*(f(xnm1mh)+f(xnm1ph));

%  Solve tridiagonal linear system.

u = A\b;

%  Compare with true solution.  Plot true and computed solution and error.
utrue = x.*(ones(n-1,1)-x);
err = norm(u-utrue)/sqrt(n),
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')