Like this? I've used: (T(1:nx,2) - T(1:nx,1))/dy = 5 and rearranged this to get T(1:nx,1). In addition I've set the y values to decrease from ny-1, instead of increasing from 2 and shifted the calculation of T(1:nx,1) to the end of the loop. Check carefully!!
for I=1:2
tic
nx = 21; %step size x -direction
ny = 21; %step size y -direction
Lx = 1;
Ly = 1;
x = linspace(0,Lx,nx);
y = linspace(0,Ly,nx);
dx = x(2)-x(1);
dy = y(2)-y(1);
k = 0;
% Intial condition
T = zeros(nx,ny);
% BC
T(1:nx,nx) = (sin(pi*x(:,1)));
%T(1:nx,1) = 5;
T(ny,1:ny) = 0;
T(1,1:ny) = (1-y(1,:));
Told = T;
tol = 1.0e-4;
error = 1;
k = k+1;
% iterative solver
iterative_solver = input('iterative_solver(1. jacobi 2. gauss) = ');
% Jacobi
if iterative_solver ==1
jacobi_iter = 1;
while(error>tol)
for i = 2:nx-1
for j = ny-1:-1:2 %%%%%%%%%%%%%%%%
T(i,j) = 0.25*(Told(i-1,j) + Told(i+1,j) + Told(i,j-1) + Told(i,j+1));
end
end
T(1:nx,1) = T(1:nx,2) - 5*dy; %%%%%%%%%%%%%
error = max(max(abs(Told-T)));
jacobi_iter = jacobi_iter+1;
Told=T;
end
contourf(x,y,T,'ShowText','on')
colorbar;
title_text=sprintf('total jacobi iteration = %d @steady state',jacobi_iter);
title(title_text)
xlabel('x-axis');
ylabel('y-axis');
end
% Gauss
if iterative_solver ==2
gs_iter =1;
while (error>tol)
for i = 2:nx-1
for j = ny-1:-1:2 %%%%%%%%%%%%%%%%%%%
T(i,j) = (0.25)*(T(i-1,j)+T(i+1,j)+T(i,j-1)+T(i,j+1));
end
end
T(1:nx,1) = T(1:nx,2) - 5*dy; %%%%%%%%%%%%%%%%%%%
error = max(max(abs(Told-T)));
gs_iter=gs_iter+1;
Told = T;
end
contourf(x,y,T,'ShowText','on')
colorbar;
xlabel('x-axis')
ylabel('y-axis')
title_text=sprintf('total gauss seidel iteration @steady state = %d',gs_iter);
title(title_text)
end
if I==1
ct_j =toc;
fprintf('%d jacobi iterations for computational time of %0.3g secondsn',k,ct_j)
elseif I==2
ct_gs = toc;
fprintf('%d gauss seidel iterations for computational time of %0.3g secondsn',k,ct_gs)
end
end
