if j == 1
if i==1
T_new(j,i,t+1) = 1/(1+2*sx+2*sy+2*sy*dy*c1)*omega*(2*sy*T_old(j+1,i,t+1)-2*sy*dy*c2+ ...
+2*sx*T_old(j,i+1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) + (1-omega)*T_old(j,i,t+1);
elseif i==i_max
T_new(j,i,t+1) = 1/(1+2*sx+2*sy+2*sy*dy*c1)*omega*(2*sy*T_old(j+1,i,t+1)-2*sy*dy*c2+ ...
+2*sx*T_old(j,i-1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) + (1-omega)*T_old(j,i,t+1);
else
T_new(j,i,t+1) = 1/(1+2*sx+2*sy+2*sy*dy*c1)*omega*(2*sy*T_old(j+1,i,t+1)-2*sy*dy*c2+ ...
sx*T_old(j,i-1,t+1)+sx*T_old(j,i+1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) + (1-omega)*T_old(j,i,t+1);
end
elseif j == j_max
if i==1
T_new(j,i,t+1) = 1/(1+2*sx+2*sy-2*sy*dy*c1)*omega*(2*sy*T_old(j-1,i,t+1)+2*sy*dy*c2+ ...
+2*sx*T_old(j,i+1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) + (1-omega)*T_old(j,i,t+1);
elseif i==i_max
T_new(j,i,t+1) = 1/(1+2*sx+2*sy-2*sy*dy*c1)*omega*(2*sy*T_old(j-1,i,t+1)+2*sy*dy*c2+ ...
+2*sx*T_old(j,i-1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) + (1-omega)*T_old(j,i,t+1);
else
T_new(j,i,t+1) = 1/(1+2*sx+2*sy-2*sy*dy*c1)*omega*(2*sy*T_old(j-1,i,t+1)+2*sy*dy*c2+ ...
sx*T_old(j,i-1,t+1)+sx*T_old(j,i+1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) + (1-omega)*T_old(j,i,t+1);
end
else
if (i==1)
T_new(j,i,t+1) = 1/(1+2*sx+2*sy)*omega*(sy*T_old(j+1,i,t+1)+sy*T_old(j-1,i,t+1)+ ...
+2*sx*T_old(j,i+1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) + (1-omega)*T_old(j,i,t+1);
elseif (i==i_max)
T_new(j,i,t+1) = 1/(1+2*sx+2*sy)*omega*(sy*T_old(j+1,i,t+1)+sy*T_old(j-1,i,t+1)+ ...
+2*sx*T_old(j,i-1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) + (1-omega)*T_old(j,i,t+1);
else
T_new(j,i,t+1) = 1/(1+2*sx+2*sy)*omega*(sy*T_old(j+1,i,t+1)+sy*T_old(j-1,i,t+1)+ ...
+sx*T_old(j,i-1,t+1)+sx*T_old(j,i+1,t+1)+T_old(j,i,t)+dt*Q(j,i,t+1)/rho/cp) + (1-omega)*T_old(j,i,t+1);
end
end
temp = temp + abs((T_new(j,i,t+1)-T_old(j,i,t+1))/T_new(j,i,t+1));
temp = temp + abs((T_new(j,i,t+1)-T_old(j,i,t+1)));
T_old(j,i,t+1) = T_new(j,i,t+1);
end
end
T_old = T_new;
change(iter) = temp/(i_max-2)/(j_max-2);
if( change(iter) < change_want )
fprintf('\n The total number of iterations excuted is %g \n',iter)
disp('Desired accuracy achieved; breaking out of main loop');
break;
end