%
% Crank-Nicolson method for the heat equation
% with homogeneous Dirichlet boundary conditions
% elementary code which do not exploit Matlab peculiarities
%
% Shows the order of the error for this method in time (wrt k, the time step)
%
clear all;
global A




mm = 10;



h = 1/mm ;
m=mm-1;
x = [h:h:1-h] ;

n = 10 ;
%interval of time [0,T]=[0,1];
T=1;
k = T/n ;


u = sin(pi*x)' ;


% u = ones(size(x))' ;
U = [0;u;0];

r=k/(2*h^2);
e=ones(m,1);
A=1/(h^2)*spdiags([e -2*e e],-1:1,m,m);

options=odeset('AbsTol',1e-10,'RelTol',1e-10);
[tt uu]=ode15s('Au',[0 1],u,options);




nexp=5;

kk=[];
for j=1:nexp
    n=n*2;
    u = sin(pi*x)' ;
    U = [0;u;0];
    k=1/n;
    kk=[kk k];

    Aleft=speye(m)-k/2*A;


% time stepping

for l=1:n,

   rhs = u+k/2*A*u ;
   u   = Aleft\rhs ;
   U   = [U,[0;u;0]] ;

end
err(j)=norm(uu(end,:)'-U(2:end-1,end));
end

figure(3)
loglog(kk, err,'-o')
hold on
loglog(kk,kk.^2,'k')
title('loglog plot off the norm 2 of the error vs the time step')
xlabel('k')
ylabel('norm of the error')
hold off
figure(1)
t = [0:k:1] ; 
mesh([0 x 1],t,U')
view(140, 30)
set(gca,'FontSize',15)
xlabel('space x')
ylabel('time t')
zlabel('u')


figure(2)
contourf([0 x 1],t,U',20)
set(gca,'FontSize',15)
colorbar;
xlabel('space x')
ylabel('time t')