clear clc
L=0.02;H=40;dx=0.01;dy=dx;dt=15
nx=uint32(L/dx+1);ny=uint32(H/dy+1);rx=dt/dx^2;ry=dt/dy^2;
k=28;alfa=12.5*10^-6;g=5*10^6;Ts=30;h=45
[X Y]=meshgrid(linspace(0,L,nx),linspace(0,H,ny));Tint=200;T=Tint*ones(ny,nx);
Tleft=0;Tright=Tint;Ttop=Tint;Tbottom=Tint;
T(:,1)=Tleft;T(:,end)=Tright;T(1,:)=Ttop;T(end,:)=Tbottom;Fo=alfa*dt*(dx^2);
time=0;n=time/dt;
err=1
tol=10^-1
while err>tol
n=n+1
Tn=T;
for i=2:nx-1
for j=2:ny-1
T(j,i)=Tn(j,i)+((Tn(j,i+1)+Tn(j,i-1)+Tn(j-1,i)+Tn(j+1,i)-4*Tn(j,i))+g*dx^2/k)*Fo
T(j,end)=Fo*((g*dx^2)/k - (2*h*(Tn(j,end) - Ts)*dx)/k - 2*Tn(j,end) + 2*Tn(j,end-1) + Tn(j,end)/Fo)
end
end
err = max(max(abs(Tn - T)))
end