%Exe 2.14
%%%%%%%%%%%

%S1:
%%

%%%%%%%% f.m %%%%%%%%%%%
function y=f(x)
y=-exp(x-1);
%%%%%%%%%%%%%%%%%%%%%

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

%S2:
%%%
n=7;h=1/(n+1);
A=spdiags(ones(n+1,1)*[-1 2 -1],[-1 0 1],n+1,n+1)
A(n+1,n)=-2
x0=linspace(0,1,n+2);x=x0(2:n+2);
b=(h^2)*f(x)'
b(n+1)=b(n+1)+2*h
v2=A\b
% Try larger:

n=70;h=1/(n+1);
A=spdiags(ones(n,1)*[-1 2 -1],[-1 0 1],n,n);
A(n,n)=1;
lam=eigs(A); sum(lam>0)==n
min(lam),max(lam)

n=70;h=1/(n+1);
A=spdiags(ones(n+1,1)*[-1 2 -1],[-1 0 1],n+1,n+1)
A(n+1,n)=-2
lam=eigs(A); sum(lam>0)==n+1
min(lam),max(lam)

%(c)
%Let's define f.m as above and exact sol u.m
%%%%%%%% f.m %%%%%%%%%%%
function y=f(x)
y=-exp(x-1);
%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%% u.m %%%%%%%%%%%%
function y=u(x)
y=exp(-1)*(exp(x)-1);

% S1
%%%

n=7;h=1/(n+1);
A=spdiags(ones(n,1)*[-1 2 -1],[-1 0 1],n,n)
A(n,n)=1
x0=linspace(0,1,n+2);x=x0(2:n+1);
b=(h^2)*f(x)'
b(n)=b(n)+h
v=A\b;v0=[0;v;v(n)+h];
clf;fplot('u',[0,1]);grid;hold on
plot(x0,v0,'or')
% Error
%%%%%%%%
Eh=max(abs(u(x0)-v0'))


%S2:
%%%
n=7;h=1/(n+1);
A=spdiags(ones(n+1,1)*[-1 2 -1],[-1 0 1],n+1,n+1)
A(n+1,n)=-2
x0=linspace(0,1,n+2);x=x0(2:n+2);
b=(h^2)*f(x)'
b(n+1)=b(n+1)+2*h
v=A\b;v0=[0;v]
clf;fplot('u',[0,1]);grid;hold on
plot(x0,v0,'or')
% Error
%%%%%%%%
Eh=max(abs(u(x0)-v0'))

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);
    A(n,n)=1;
    x0=linspace(0,1,n+2);x=x0(2:n+1);
    b=(h^2)*f(x)';
    b(n)=b(n)+h;
    v=A\b;v0=[0;v;v(n)+h];
    EH=[EH,max(abs(u(x0)-v0'))];
    clf;fplot('u',[0,1]);grid;hold on
    plot(x0,v0,'or'); pause
end;

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]']

N=80;EH=[];
for n=[5 10 20 40 80];
    h=1/(n+1);
    A=spdiags(ones(n+1,1)*[-1 2 -1],[-1 0 1],n+1,n+1);
    A(n+1,n)=-2;
    x0=linspace(0,1,n+2);x=x0(2:n+2);
    b=(h^2)*f(x)';
    b(n+1)=b(n+1)+2*h;
    v=A\b;v0=[0;v];
    EH=[EH,max(abs(u(x0)-v0'))];
    clf;fplot('u',[0,1]);grid;hold on
    plot(x0,v0,'or'); pause
end;

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    ratio']
[N' H' EH' [1,ratio]']


% Exe 2.17
%%%%%%%%%%%

%%%%%%%%%%% f.m %%%%%%%%%%
function y=f(x)
y=100*exp(-10*x);
%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%% u.m %%%%%%%%%%
function y=u(x)
y=1+x*exp(-10)-x-exp(-10*x);
%%%%%%%%%%%%%%%%%%%%%%%%%

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'))];
    clf;fplot('u',[0,1]);grid;hold on
    plot(x0,v0,'or'); pause
end;

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]']

errest=(H.^2)*max(abs(10000*exp(-10*x)))/96
[N' H' EH' errest']


>> [N' H' EH' errest']

ans =

    5.0000    0.1667    0.1305    2.5575
   10.0000    0.0909    0.0438    0.7609
   20.0000    0.0476    0.0125    0.2088
   40.0000    0.0244    0.0033    0.0548
   80.0000    0.0123    0.0008    0.0140


