Matlab program with the Crank-Nicholson method

Matlab program with the Crank-Nicholson method for the diffusion equation

https://www.youtube.com/watch?v=HiCdEG7uw7E

%Heat diffusion in one dimension wire within

%Parameters to define the heat equation
cond = 1/2;   % Heat conductivity parameter
L = 1.;       % Length of the wire
T = 1.;       % Final Time

%Parametes for Crank-Nocholson
Nt = 2500;          % Number of time steps
Dt = T/Nt;
Nx = 50;            % Number of space steps
Dx = L/Nx;
b = cond*Dt/(Dx*Dx);% Alpha Value

%Initial Condition
for n = 1:Nx-1
    x(n) = (n-1)*Dx;
    u(n,1) = sin(pi*x(n));
end

%Boundary Condition
for j = 1:Nt+1
    u(1,j) = 0.;  % Boundary is 0 for any moment
    u(Nx+1,j) = 0.;
    t(j) = (j-1)*Dt;
end

%Defining the M_right and M_left matrices in Crank method
% As mmr 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;                      %Abpve
MMl = diag(bbl,0) + diag(aal,-1)+diag(ccl,1);


aar(1:Nx-2) = b;                      %Below the main diagonal in MML
bbr(1:Nx-1) = 2.-2.*b;                 % The main diagonal in MML
ccr(1:Nx-2) = b;                      %Abpve
MMR = diag(bbr,0) + diag(aar,-1)+diag(ccr,1);

%Crank
for j = 1:Nt
    u(2:Nx,j+1) = inv(MML) * MMR * u(2:Nx,j);
end

figure(1);
plot(x,u(:,1),'-b');
plot(x,u(:,round(Nt/100)),'-g');
plot(x,u(:,round(Nt/10)),':b');
plot(x,u(:,Nt),'-.r');

xlabel('X');
ylabel('T');

figure(2);
mesh(t,x,u);
xlabel('t');
ylabel('x');