Two point boundary value problems using FEM

[EEHJ] 

ff.m   % There is a bug in quad, file name can't be f.m
---
function y=ff(x)
y=(3*x+x.^2).*exp(x);
----

%Test:
%%%%%%%

m=8;
x=linspace(0,1,m+1);M=length(x);
h=x(2:M)-x(1:M-1)
inind=2:M-1;   % Inner node indices
xin=x(inind);%%hin=h(inind);
d=1./h(1:M-2)+1./h(2:M-1);
ud=-1./h(2:M-2);ld=-1./h(1:M-3);
A=spdiags([[ud,0]', d', [0,ld]'],[-1 0 1],M-2,M-2)
full(A)


---------
function y=xxf(x)
y=x.*ff(x);
---------
clear b
...


U=A\b

U=[0;U;0];
plot(x,U,'or')


% 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];

% Exact sol (from Maple session)
% u.m
function y=u(x)
y=x.*exp(x).*(1-x);
%%%%%%%%%%%%%
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)*f(x)';
    v=A\b;v0=[0;v;0];
    EH=[EH,max(abs(u(x0)-v0'))];
end;

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

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

ans =

         n                  h                 Eh
   5.00000000000000   0.16666666666667   0.00588534055622
  10.00000000000000   0.09090909090909   0.00178472554444
  20.00000000000000   0.04761904761905   0.00049104625305
  40.00000000000000   0.02439024390244   0.00012883077371
  80.00000000000000   0.01234567901235   0.00003301624304


nn=length(EH); ratio=EH(2:nn)./EH(1:nn-1)

[N' H' EH' [1,ratio]']

>> nn=length(EH); ratio=EH(2:nn)./EH(1:nn-1)

ratio =

   0.30324932387437   0.27513824440864   0.26235975309927   0.25627605958246

>> [N' H' EH' [1,ratio]']                   

ans =
         n                  h                  Eh                 ratio
                                                            (first is for fill)
   5.00000000000000   0.16666666666667   0.00588534055622   1.00000000000000 
  10.00000000000000   0.09090909090909   0.00178472554444   0.30324932387437
  20.00000000000000   0.04761904761905   0.00049104625305   0.27513824440864
  40.00000000000000   0.02439024390244   0.00012883077371   0.26235975309927
  80.00000000000000   0.01234567901235   0.00003301624304   0.25627605958246