% Here's the "oldfashioned" FD-solution for comparison:
%%%%%%%

n=7;h=1/(n+1);
A=spdiags(ones(n,1)*[-1 2 -1],[-1 0 1],n,n)
x0=linspace(0,1,n+2);x=x0(2:n+1);
b=(h^2)*ff(x)';%% b0=[ff(0);b;ff(1)];

%Solve:
%%%%%%%

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

%%%%%%%%%%%%%
clf;fplot('u',[0,1]);grid;hold on
plot(x0,v0,'o')

%%axis([-.1 1.1 -.1 .45])

% Error
%%%%%%%%
Eh=max(abs(u(x0)-v0'))

clear EH
%%N=5*(1:5);H=1./(N+1)
N=80;EH=[];
for n=[5 10 20 40 80];
    h=1/(n+1);
    A=spdiags(ones(n,1)*[-1 2 -1],[-1 0 1],n,n);
    x0=linspace(0,1,n+2);x=x0(2:n+1);
    b=(h^2)*ff(x)';
    v=A\b;v0=[0;v;0];
    EH=[EH,max(abs(u(x0)-v0'))];
end;

format long;
N=[5 10 20 40 80];H=1./(N+1);
[N' H' EH']

nn=length(EH); ratio=EH(2:nn)./EH(1:nn-1)
[N' H' EH' [1,ratio]']

hold off;
