%  Heat diffusion in one dimensional wire within the Crank-Nicholson method
clear;

% Parameters to define the heat equation and the range in space and time
cond = 1/2;		%  Heat conductivity parameter
L = 1.;  		%  Length of the wire
T =1.;			%  Final time

%  Parameters needded to solve the equation within the  Crank-Nicholson method
Nt = 2500;					%  Number of time steps
Dt = T/Nt;                  %  Time step
Nx = 50;					%  Number of space steps
Dx = L/Nx;                  %  Space step
b = cond*Dt/(Dx*Dx);		%  Alpha parameter in the finite-difference implementation



%  The initial condition: the initial temperature of the wire, in this example a sin function
for n = 1:Nx+1
        x(n) =(n-1)*Dx;     % We also define vector x, due to the space discretization
        u(n,1) =sin(pi*x(n));        
end
          
%  Boundary conditions: the temperature of the wire at the boundaries at any time (in this example T=0)
for j=1:Nt+1
        u(1,j) = 0.;
        u(Nx+1,j) = 0.;
        t(j) = (j-1)*Dt; % We also define vector t, due to the time discretization
end


% Defining the  M_right (MMr) and M_left (MMl) matrices in the Crank-Nicholson method
% As MMr and MMl matrices are tridiagonal we first define the three diagonals and then build the matrices

aal(1:Nx-2)=-b;                             % Below the main diagonal in MMl
bbl(1:Nx-1)=2.+2.*b;                        % The main diagonal in MMl
ccl(1:Nx-2)=-b;                             % Above the main diagonal in MMl
MMl=diag(bbl,0)+diag(aal,-1)+diag(ccl,1);   % Building the MMl matrix

aar(1:Nx-2)=b;                              % Below the main diagonal in MMr
bbr(1:Nx-1)=2.-2.*b;                        % The main diagonal in MMr
ccr(1:Nx-2)=b;                              % Above the main diagonal in MMr
MMr=diag(bbr,0)+diag(aar,-1)+diag(ccr,1);   % Building the MMr matrix


% Implementation of the Crank-Nicholson method
for j=1:Nt                                  % Time Loop   
   u(2:Nx,j+1)=inv(MMl)*MMr*u(2:Nx,j);      % Implementation of the Crank-Nicholson method to get u at time j
end

% Graphical representation of the temperature at different selected times
figure(1)
plot(x,u(:,1),'-b',x,u(:,round(Nt/100)),'--g',x,u(:,round(Nt/10)),':b',x,u(:,Nt),'-.r')
title('Temperature of the wire within the Crank-Nicholson method')
xlabel('X')
ylabel('T')

% 3D graphical representation of the temperature as a function of time
figure(2)
mesh(t,x,u)

title('Temperature of the wire within the Crank-Nicholson method')
xlabel('t')
ylabel('x')
