%% EXAMPLE 3: More exhaustive testing

exact = @(x) sin(4*pi*x);
u0=0; u1=0;
a = @(x)   x.^2+1;
b = @(x)   1./(1+x);
c = @(x)   3*x;
f = @(x)   -(x.^2+1).*(-16*pi^2*sin(4*pi*x))+...
            4*pi*cos(4*pi*x)./(1+x)+...
            3*x.*sin(4*pi*x);
M=2.^(2:12);    % [4 8 16 32 ....]
for i=1:length(M)
    [u,x]=fdBVP1d(a,b,c,f,u0,u1,M(i));
    error(i)=max(abs(u-exact(x)));
end
disp(error)
h=1./M;
plot(-log10(h),-log10(error),'o-'), xlabel('-log_{10}(h)'), ylabel('-log_{10}(Error_h)')
ecr=log10(error(1:end-1)./error(2:end))./log10(h(1:end-1)./h(2:end));
       % ecr= estimated convergence rate
disp(ecr)

return

%% EXAMPLE 3: Oscillating exact solution

M=input('Give M: ');
exact = @(x) sin(4*pi*x);
u0=0; u1=0;
a = @(x)   x.^2+1;
b = @(x)   1./(1+x);
c = @(x)   3*x;
f = @(x)   -(x.^2+1).*(-16*pi^2*sin(4*pi*x))+...
            4*pi*cos(4*pi*x)./(1+x)+...
            3*x.*sin(4*pi*x);
[u,x]=fdBVP1d(a,b,c,f,u0,u1,M);
plot(x,u,'-b'), hold on
plot(x,exact(x),'r'), hold off
max(abs(u-exact(x)))

return


%% EXAMPLE 2: Exact solution is polynomial of degree 2
% It has to be computed exactly

exact = @(x) 1+x.^2+3*x;
u0=1; u1=5;
a = @(x)   x.^2+1;
b = @(x)   1./(1+x);
c = @(x)   3*x;
f = @(x)   -(x.^2+1)*2+(2*x+3)./(1+x)+3*x.*(1+x.^2+3*x);
[u,x]=fdBVP1d(a,b,c,f,u0,u1,100);
plot(x,u,'-'), hold on
plot(x,exact(x),'r'), hold off
max(abs(u-exact(x)))

return


%% EXAMPLE 1: Straight line solution

a = @(x)   ones(size(x));    % 1+0*x also works
b = @(x)   0*x;
c = @(x)   0*x;
f = @(x)   0*x;
u0=1; u1=2;
exact = @(x)  u0+(u1-u0)*x;
[u,x]=fdBVP1d(a,b,c,f,u0,u1,10);
plot(x,u,'o-')     % You should add tags to the axes
max(abs(u-exact(x)))

return