%  Inline functions define the conductivity (kappa), the source
%  function (f), and the true solution (u).

kappa = inline('1+x^2');
f = inline('2*(3*x^2-2*x+1)');
u = inline('(1-x)^2');

%  Ask user for problem size.

n = input(' Enter number of subintervals n: ');
h = 1/n;
A = sparse(n,n);  % Store A as a sparse matrix for efficiency.
F = zeros(n,1);

%  Set rows 1 through n-1 of A and F.

for i=1:n-1,
  xi = i*h;
  kappaimh = kappa(xi-h/2);  kappaiph = kappa(xi+h/2);
  A(i,i) = -(kappaimh+kappaiph);
  A(i,i+1) = kappaiph;
  if i > 1, A(i,i-1) = kappaimh; end;
  F(i) = f(xi);
  if i==1, F(i) = F(i) - kappaimh/h^2; end;
end;

%  Set row n of A and F.

xn = n*h;
kappanmh = kappa(xn-h/2);  kappanph = kappa(xn+h/2);
A(n,n) = -(kappanmh+kappanph);
A(n,n-1) = kappanmh+kappanph;
F(n) = f(xn);

%  Remember to multiply A by h^{-2}.

A = A/h^2;

%  Solve linear system.

uapprox = A\F;

%  Compare to true solution.

utrue = zeros(n,1);
for i=1:n, xi = i*h; utrue(i) = u(xi); end;
err2 = sqrt(h)*norm(utrue-uapprox),
errinf = max(abs(utrue-uapprox)),

%  Plot error vector.
plot([0:h:1]',[0;utrue-uapprox])
xlabel('x'), ylabel('error')