f = inline('2*(3*x^2 - 2*x + 1)');  % Set f so that soln is (1-x)^2.

n = input(' Enter n: ');  % Ask user for number of subintervals.
h = 1/n,
x = [h:h:1]';  % Set nodes.

%  Initialize array A and rhs vector b.
A = zeros(n,n); b = zeros(n,1);

%  Set matrix and right-hand side.
for i=1:n,
  A(i,i) = -(1/h^2)*((x(i)+h/2)^2 + (x(i)-h/2)^2 + 2);
  if i < n, A(i,i+1) = (1/h^2)*(1 + (x(i)+h/2)^2); end;
  if i > 1, A(i,i-1) = (1/h^2)*(1 + (x(i)-h/2)^2); end;
  if i==n, A(i,i-1) = -A(i,i); end;
  b(i) = f(x(i));
  if i==1, b(i) = b(i) - (1/h^2)*(1 + (x(i)-h/2)^2); end;
end;

%  Solve for u.
u = A\b;

%  Set true solution utrue and compare to u.
utrue = (1 - x).^2;
errL2 = norm(utrue-u)/sqrt(n),