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 %d\n', iter);
    end
end

Iteration 100
Iteration 200
Iteration 300
Iteration 400
Iteration 500
Iteration 600
Iteration 700
Iteration 800
Iteration 900
Iteration 1000

% 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');

[Torsten]

Where do you update Psi in your code ? I guess it's the solution of Laplace's equation ?

[Ruthra]

nabla^{2} \Psi=-\zeta