A convergence prcedure that doesn't work

Hello!!

I use the Crank-Nicolson numerical method to calculate temperatures on a metalic plate(aluminium) 10cmx30cm.

I assume that the left vertical side is always in 0 degrees Celsius and the right vertical side in 100 degrees.Both horizontal sides are insulated. So we have a heat diffusion problem. Crank-Nicolson numerical method, says in few words, we divide the plate in intervals using points, each point 's temperature is calculated by using the temperatures of its adjacent points in x,y directions.

T(i,j) = 1/2*a*dt*[E(i-1,j)-2E(i,j)+E(i+1,j)+ T(i-1,j)-2T(i,j)+T(i+1,j)]/dx^2 +

+ 1/2*a*dt*[E(i,j-1)-2E(i,j)+E(i,j+1)+ T(i,j-1)-2T(i,j)+T(i,j+1)]/dy^2 +

+E(i,j)

E is the temperature for the previous time point(we have also divided time in dt intervals) For t=0 E is known(20 degrees)

For the uknown values of T we use an arbitary value(lets say 1) and by repeating the calculation for all the points' temperatures T values must converge on the right ones.

I wrote a very simple algorithm that calculates T(i,j) for the first second(dt=1) when the plate is divided by 2cm intervals (dx=dy=2cm this means i have 6 points on y and 16 on x direction) All i care about is the convergence procedure. It works fine!!

THE PROBLEM IS, it doesn't work if i just use dx=dy=1cm which means a 11X31 matrix Here is my code that works

clear all close all clc

E=zeros(6,16); T=ones(6,16);

dx=2; dy=2; dt=1; a=0.8418; for i=1:6 for j=2:15 E(i,j)=20; end end

for i=1:6 E(i,1)=0; E(i,16)=100; T(i,1)=0; T(i,16)=100; end

for m=1:50 for j=2:15 T(1,j)=(4/3)*T(2,j)-(1/3)*T(3,j); T(6,j)=(4/3)*T(5,j)-(1/3)*T(4,j); end for i=2:5 for j=2:15 T(i,j)=(1/2)*dt*a*(E(i-1,j)-2*E(i,j)+E(i+1,j)+T(i-1,j)-2*T(i,j)+T(i+1,j))/dx^2+(1/2)*dt*a*(E(i,j-1)-2*E(i,j)+E(i,j+1)+T(i,j-1)-2*T(i,j)+T(i,j+1))/dy^2 + E(i,j);

    end
end
end
for j=2:15
        T(1,j)=(4/3)*T(2,j)-(1/3)*T(3,j); 
        T(6,j)=(4/3)*T(5,j)-(1/3)*T(4,j); 
end

and here is the same code that doesn't work

clear all close all clc

E=zeros(11,31); T=ones(11,31);

dx=1; dy=1; dt=1; a=0.8418; for i=1:11 for j=2:30 E(i,j)=20; end end

for i=1:11 E(i,1)=0; E(i,31)=100; T(i,1)=0; T(i,31)=100; end

for m=1:50 for j=2:30 T(1,j)=(4/3)*T(2,j)-(1/3)*T(3,j); T(11,j)=(4/3)*T(10,j)-(1/3)*T(9,j); end for i=2:10 for j=2:30 T(i,j)=(1/2)*dt*a*(E(i-1,j)-2*E(i,j)+E(i+1,j)+T(i-1,j)-2*T(i,j)+T(i+1,j))/dx^2 + (1/2)*dt*a*(E(i,j-1)-2*E(i,j)+E(i,j+1)+T(i,j-1)-2*T(i,j)+T(i,j+1))/dy^2 + E(i,j);

    end
end
end
for j=2:30
        T(1,j)=(4/3)*T(2,j)-(1/3)*T(3,j);
        T(11,j)=(4/3)*T(10,j)-(1/3)*T(9,j);
end