% The FEM solution
npts = 7;
x=linspace(0,1,npts);
h=x(2:npts)-x(1:npts-1);

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

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

ksi = A \ b';

% Here's the "oldfashioned" FD-solution for comparison:
%%%%%%%
npts = 5;
h=1/(npts+1);
A=spdiags(ones(npts, 1)*[-1 2 -1], [-1 0 1], npts, npts);
x0=linspace(0, 1, npts+2); x=x0(2:npts+1);
b=(h^2)*ff(x)'; %% b0=[f(0);b;f(1)];

%Solve:
%%%%%%%

fd_v=A\b; % v0=[0;v;0];

% Exact sol (from Maple session)
u = x.*exp(x).*(1-x);

[x; u; ksi'; fd_v']