halfh = [0.2012 0.0673 0.0493 0.0417 0.0376 0.0353 0.0341 0.0335];
h = [halfh halfh(length(halfh):-1:1)];
npts = length(h) + 1;

x = 0;
for k=1:length(h)
  x(k+1) = x(k) + h(k);
end

d=1./h(1:npts-2) + 1./h(2:npts-1);
ud=-1./h(2:npts-2); ld=-1./h(1:npts-3);
A = spdiags([[ud, 0]', d', [0, ud]'], [-1 0 1], npts-2, npts-2);

clear b
for i=1:npts-2;
  b(i)=(quad8('xf2', x(i), x(i+1))-x(i)*quad8('f2', x(i), x(i+1)))/h(i);
  b(i)=b(i)+(x(i+2)*quad8('f2', x(i+1), x(i+2)) - ...
  quad8('xf2', x(i+1), x(i+2)))/h(i+1);
end

ksi = [0; A \ b'; 0];

% Compute the exact solution
xx = 0:0.02:1;
ts = exp(-5*xx.*xx + 5*xx - 5/4) - exp(-5/4);

plot(xx, ts, x, ksi, 'o');

