As I already wrote: We don't like to conclude from the code what you try to do. So please clarify:
What is your solution domain ? What is the equation you are trying to solve (maybe different for different parts of the domain ?) ? What are your boundary conditions ?
Look at the below code. I introduced z to count the number of nodes that are addressed in your loop. These are 40. But your grid has 42 nodes - thus two of them are never changed and their value of 1 persists.
%%Parámetros geométricos del dominio
W=2; %Grosor (metros)
A=3; %Anchura (en metros)
H=2.5; %Altura (en metros)
Nx= 7; %Número de nodos en la dirección x
Ny=6; %Número de nodos en la dirección y
dx=A/(Nx-1); %Distancia de elemento finito en dirección x (metros)
dy=H/(Ny-1); %Distancia de elemento finito en direccion y (metros)
%%Condiciones iniciales y de frontera
To=1000; %Temperatura interna del horno (en C)
T = ones(Nx,Ny);
Tinf=20; %Temperatura externa (en C)
ho=100; %Coeficiente de transferencia convectiva interior (en W/m^2K)
hinf=10; %Coeficiente de transferencia convectiva exterior (en W/m^2K)
q_dot=8000; %Calor generado por unidad volumétrica
K=1; %Conductividad del material (en W/mK)
epsilon=1e-4;
Error=5;
max_iter=15000;
%%Computación
iter=0;
Told=T;
z=0;
%while (Error > epsilon && iter<max_iter)
iter = iter + 1;
%disp(iter);
for j = 1:Ny
for i = 1:Nx
if i >= 2 && i <= Nx-1 %Pared superior (hinf y Tinf)
z=z+1;
if j == Ny
T(i, j) = (T(i - 1, j) + T(i + 1, j) + 2 * T(i, j - 1) + (q_dot * dx^2) / K + (2 * hinf * dx * Tinf) / K) / (4 + 2 * hinf * dx / K);
elseif j >= 2 && j <= Ny-1 %Nodos internos (inner nodes)
T(i, j) = (T(i - 1, j) + T(i + 1, j) + T(i, j - 1) + T(i, j + 1) - 4 * T(i, j)) * K / q_dot + T(i, j);
end
elseif i>round(Nx*(4/7)) && i<Nx %Pared inferior (bottom)
z=z+1;
if j==1
T(i,j)=(2*T(i,j+1)+T(i-1,j)+T(i+1,j)+(q_dot*dx^2)/K)/4;
end
elseif i == round(Nx * (4 / 7)) && j>=1 && j<= round(Ny*0.5) %Pared vertical interior (ho y To)
z=z+1;
if j == 1
T(i, j) = (T(i, j + 1) + T(i + 1, j) + (q_dot * dx^2) / (2 * K) + (ho * dx * To) / K)/(2+ho*dx/K);
elseif j == round(Ny * 0.5)
T(i, j) = (T(i, j + 1) + 2 * T(i + 1, j) + T(i, j - 1) + T(i - 1, j) + (3 * q_dot * dx^2) / (2 * K) + (2 * ho * dx * To) / K) / (6 + (2 * ho * dx) / K);
else
T(i, j) = ((T(i, j - 1) + T(i + 1, j) + 2 * T(i + 1, j) + (q_dot * dx^2) / K) + (2 * ho * dx * To) / K) / (4 + 2 * ho * dx / K);
end
elseif i > 1 && i < round(Nx * 4 / 7) %Pared horizontal interior (ho y T o)
z=z+1;
if j == Ny * 0.5
T(i, j) = ((T(i - 1, j) + T(i + 1, j) + 2 * T(i, j + 1) + (q_dot * dx^2) / K) + (2 * ho * dx * To) / K) / (4 + 2 * ho * dx / K);
end
elseif i == Nx && 1<=j && Ny>=j %Pared derecha
z=z+1;
if j == Ny
T(i, j) = (T(i, j - 1) + T(i - 1, j) + (q_dot * dx^2) / (2 * K) + (hinf * dx * Tinf) / K)/(2+hinf*dx/K);
elseif j == 1
T(i, j) = ((q_dot * dx^2) / (2 * K) + T(i - 1, j) + T(i, j + 1));
else
T(i, j) = (2 * T(i - 1, j) + T(i, j - 1) + T(i, j + 1) + (q_dot * dx^2) / K) / 4;
end
elseif i == 1 && j>=round(0.5*Ny) && j<=Ny %Pared Izquierda
z = z+1;
if j == round(Ny * 0.5)
T(i, j) = (T(i, j + 1) + T(i + 1, j) + (q_dot * dx^2) / (2 * K) + (ho * dx * To) / K)/(2+ho*dx/K);
elseif j == Ny
T(i, j) = (T(i, j - 1) + T(i + 1, j) + (q_dot * dx^2) / (2 * K) + (hinf * dx * To) / K)/(2+hinf*dx/K);
else
T(i, j) = (2 * T(i + 1, j) + T(i, j - 1) + T(i, j + 1) + (q_dot * dx^2) / K) / 4;
end
end
end
end
Error = max(max(abs(T-Told)));
%disp(Error);
T(1:3,1)=To;
T(1:3,2)=To;
Told=T;
z
%end
disp(T);
