anyone correct the below code for solving energy, vorticity eqn using the alternative direct implicit fdm method for which i should get the values of velocitiy at each grid i

% MATLAB code to solve the coupled temperature and vorticity equations using ADI FDM
% Parameters
Pr = 0.733; % Prandtl number
A = 1; % Length of the domain in X direction
Ha = 10; % Hartmann number (can change this to study its influence)
Gr = 2e4; % Grashof number
Nx = 41; % Number of grid points in X direction
Ny = 41; % Number of grid points in Y direction
dx = A/(Nx-1); % Grid spacing in X
dy = 1/(Ny-1); % Grid spacing in Y
dt = 0.001; % Time step size
max_iter = 1000; % Maximum iterations for time-stepping
% Initialize arrays
T = rand(Nx, Ny); % Temperature field
zeta = rand(Nx, Ny); % Vorticity field
Psi = rand(Nx, Ny); % Stream function
U = rand(Nx, Ny); % X-component of velocity
V = rand(Nx, Ny); % Y-component of velocity
% Boundary conditions for Psi and T
T(:,1) = -1; % At Y = 0, T = -1
T(:,end) = 1; % At Y = 1, T = 1
Psi(:,1) = 0; Psi(:,end) = 0; % Stream function at Y boundaries
Psi(1,:) = 0; Psi(end,:) = 0; % Stream function at X boundaries
% Time-stepping loop
for iter = 1:max_iter
% Compute velocities U and V from stream function Psi
for i = 2:Nx-1
for j = 2:Ny-1
U(i,j) = (Psi(i,j+1) – Psi(i,j-1))/(2*dy); % U = dPsi/dY
V(i,j) = -(Psi(i+1,j) – Psi(i-1,j))/(2*dx); % V = -dPsi/dX
end
end

% Solve for temperature T using ADI method (X-direction first)
for j = 2:Ny-1
for i = 2:Nx-1
T(i,j) = T(i,j) + dt * ( …
-(U(i,j) * (T(i+1,j) – T(i-1,j))/(2*dx)) – … % Advection term (X)
(V(i,j) * (T(i,j+1) – T(i,j-1))/(2*dy)) + … % Advection term (Y)
(1/Pr) * ((T(i+1,j) – 2*T(i,j) + T(i-1,j))/(dx^2) + … % Diffusion (X)
(T(i,j+1) – 2*T(i,j) + T(i,j-1))/(dy^2)) ); % Diffusion (Y)
end
end

% Solve for vorticity zeta using ADI method (Y-direction next)
for j = 2:Ny-1
for i = 2:Nx-1
zeta(i,j) = zeta(i,j) + dt * ( …
-(U(i,j) * (zeta(i+1,j) – zeta(i-1,j))/(2*dx)) – … % Advection (X)
(V(i,j) * (zeta(i,j+1) – zeta(i,j-1))/(2*dy)) + … % Advection (Y)
(Gr/2) * (T(i,j+1) – T(i,j))/(dy) + … % Buoyancy effect
((zeta(i+1,j) – 2*zeta(i,j) + zeta(i-1,j))/(dx^2) + … % Diffusion (X)
(zeta(i,j+1) – 2*zeta(i,j) + zeta(i,j-1))/(dy^2)) – … % Diffusion (Y)
Ha^2 * (V(i+1,j) – V(i-1,j))/(2*dx) ); % Hartmann term
end
end

% Enforce boundary conditions on Psi
Psi(:,1) = 0; Psi(:,end) = 0; % At Y boundaries
Psi(1,:) = 0; Psi(end,:) = 0; % At X boundaries

% Enforce boundary conditions on T
T(:,1) = -1; % At Y = 0
T(:,end) = 1; % At Y = 1
T(1,:) = T(2,:); % Neumann BC at X = 0
T(end,:) = T(end-1,:); % Neumann BC at X = A

% Optional: Display or store the solution at certain intervals
if mod(iter, 100) == 0
fprintf(‘Iteration %dn’, iter);
end
end
% Final solution display – Single graph with subplots for T and Psi
figure;
subplot(1,2,1);
contourf(linspace(0,A,Nx), linspace(0,1,Ny), T’, 20); colorbar;
title(‘Temperature Distribution T’);
xlabel(‘X’); ylabel(‘Y’);
subplot(1,2,2);
streamslice(linspace(0,A,Nx), linspace(0,1,Ny), U’, V’); % Plot streamlines for velocity
title(‘Streamlines – Fluid Flow’);
xlabel(‘X’); ylabel(‘Y’);% MATLAB code to solve the coupled temperature and vorticity equations using ADI FDM
% Parameters
Pr = 0.733; % Prandtl number
A = 1; % Length of the domain in X direction
Ha = 10; % Hartmann number (can change this to study its influence)
Gr = 2e4; % Grashof number
Nx = 41; % Number of grid points in X direction
Ny = 41; % Number of grid points in Y direction
dx = A/(Nx-1); % Grid spacing in X
dy = 1/(Ny-1); % Grid spacing in Y
dt = 0.001; % Time step size
max_iter = 1000; % Maximum iterations for time-stepping
% Initialize arrays
T = rand(Nx, Ny); % Temperature field
zeta = rand(Nx, Ny); % Vorticity field
Psi = rand(Nx, Ny); % Stream function
U = rand(Nx, Ny); % X-component of velocity
V = rand(Nx, Ny); % Y-component of velocity
% Boundary conditions for Psi and T
T(:,1) = -1; % At Y = 0, T = -1
T(:,end) = 1; % At Y = 1, T = 1
Psi(:,1) = 0; Psi(:,end) = 0; % Stream function at Y boundaries
Psi(1,:) = 0; Psi(end,:) = 0; % Stream function at X boundaries
% Time-stepping loop
for iter = 1:max_iter
% Compute velocities U and V from stream function Psi
for i = 2:Nx-1
for j = 2:Ny-1
U(i,j) = (Psi(i,j+1) – Psi(i,j-1))/(2*dy); % U = dPsi/dY
V(i,j) = -(Psi(i+1,j) – Psi(i-1,j))/(2*dx); % V = -dPsi/dX
end
end

% Enforce boundary conditions on Psi
Psi(:,1) = 0; Psi(:,end) = 0; % At Y boundaries
Psi(1,:) = 0; Psi(end,:) = 0; % At X boundaries

% Enforce boundary conditions on T
T(:,1) = -1; % At Y = 0
T(:,end) = 1; % At Y = 1
T(1,:) = T(2,:); % Neumann BC at X = 0
T(end,:) = T(end-1,:); % Neumann BC at X = A

% Optional: Display or store the solution at certain intervals
if mod(iter, 100) == 0
fprintf(‘Iteration %dn’, iter);
end
end
% Final solution display – Single graph with subplots for T and Psi
figure;
subplot(1,2,1);
contourf(linspace(0,A,Nx), linspace(0,1,Ny), T’, 20); colorbar;
title(‘Temperature Distribution T’);
xlabel(‘X’); ylabel(‘Y’);
subplot(1,2,2);
streamslice(linspace(0,A,Nx), linspace(0,1,Ny), U’, V’); % Plot streamlines for velocity
title(‘Streamlines – Fluid Flow’);
xlabel(‘X’); ylabel(‘Y’); % MATLAB code to solve the coupled temperature and vorticity equations using ADI FDM
% Parameters
Pr = 0.733; % Prandtl number
A = 1; % Length of the domain in X direction
Ha = 10; % Hartmann number (can change this to study its influence)
Gr = 2e4; % Grashof number
Nx = 41; % Number of grid points in X direction
Ny = 41; % Number of grid points in Y direction
dx = A/(Nx-1); % Grid spacing in X
dy = 1/(Ny-1); % Grid spacing in Y
dt = 0.001; % Time step size
max_iter = 1000; % Maximum iterations for time-stepping
% Initialize arrays
T = rand(Nx, Ny); % Temperature field
zeta = rand(Nx, Ny); % Vorticity field
Psi = rand(Nx, Ny); % Stream function
U = rand(Nx, Ny); % X-component of velocity
V = rand(Nx, Ny); % Y-component of velocity
% Boundary conditions for Psi and T
T(:,1) = -1; % At Y = 0, T = -1
T(:,end) = 1; % At Y = 1, T = 1
Psi(:,1) = 0; Psi(:,end) = 0; % Stream function at Y boundaries
Psi(1,:) = 0; Psi(end,:) = 0; % Stream function at X boundaries
% Time-stepping loop
for iter = 1:max_iter
% Compute velocities U and V from stream function Psi
for i = 2:Nx-1
for j = 2:Ny-1
U(i,j) = (Psi(i,j+1) – Psi(i,j-1))/(2*dy); % U = dPsi/dY
V(i,j) = -(Psi(i+1,j) – Psi(i-1,j))/(2*dx); % V = -dPsi/dX
end
end

% Enforce boundary conditions on Psi
Psi(:,1) = 0; Psi(:,end) = 0; % At Y boundaries
Psi(1,:) = 0; Psi(end,:) = 0; % At X boundaries

% Enforce boundary conditions on T
T(:,1) = -1; % At Y = 0
T(:,end) = 1; % At Y = 1
T(1,:) = T(2,:); % Neumann BC at X = 0
T(end,:) = T(end-1,:); % Neumann BC at X = A

Cart

Cart