solution for steady flow at Re = 10. in coursera i took mathematics for engineers.for that i need to solve an assignment. For this code i cant get correct stream & vorticity

flow_around_cylinder_steady

function [psi, omega] = flow_around_cylinder_steady
Re=10;

%%%%% define the grid %%%%%
n=101; m=101; % number of grid points
N=n-1; M=m-1; % number of grid intervals
h=pi/M; % grid spacing based on theta variable
xi=(0:N)*h; theta=(0:M)*h; % xi and theta variables on the grid

%%%%% Initialize the flow fields %%%%%
psi=zeros(n,m);
omega=zeros(n,m);
psi(n,:)=0 % Free stream boundary condition (psi = 0 at the top edge)

%%%%% Set relax params, tol, extra variables %%%%%
r_psi=1.8; % Set the relaxation parameter here, psi equation
r_omega=0.9; % Set the relaxation parameter here, omega equation
delta=1.e-08; % error tolerance
error=2*delta; % initialize error variable

%%%%% Add any additional variable definitions here %%%%%
…
…
%%%%% Main SOR Loop %%%%%
while (error > delta)
psi_old = psi; omega_old = omega;
for i=2:n-1
for j=2:m-1
psi(i,j)=… % (1 – r_psi) * psi_old(i, j) + …
r_psi * 0.5 * (psi_old(i+1, j) + psi_old(i-1, j) + …
psi_old(i, j+1) + psi_old(i, j-1) – h^2 * omega(i, j));
end
end
error_psi=max(abs(psi(:)-psi_old(:)));
omega(1,:)=0 % Boundary condition at theta = 0
for i=2:n-1
for j=2:m-1
omega(i,j)=… % (1 – r_omega) * omega_old(i, j) + …
r_omega * ( (omega_old(i+1, j) + omega_old(i-1, j) + …
omega_old(i, j+1) + omega_old(i, j-1) – …
(4 – (2/h^2)) * omega_old(i, j)) / (4 – (2/h^2)) );
end
end
error_omega=max(abs(omega(:)-omega_old(:)));
error=max(error_psi, error_omega);
end
plot_Re10(psi);
end
Please refer to page 59 (64/69) of the PDF document for the ‘plot_Re10’ custom plotting function.
https://www.math.hkust.edu.hk/~machas/flow-around-a-cylinder.pdfflow_around_cylinder_steady

function [psi, omega] = flow_around_cylinder_steady
Re=10;

%%%%% Initialize the flow fields %%%%%
psi=zeros(n,m);
omega=zeros(n,m);
psi(n,:)=0 % Free stream boundary condition (psi = 0 at the top edge)

function [psi, omega] = flow_around_cylinder_steady
Re=10;

%%%%% Initialize the flow fields %%%%%
psi=zeros(n,m);
omega=zeros(n,m);
psi(n,:)=0 % Free stream boundary condition (psi = 0 at the top edge)

Cart

Cart