Why is the Omega result not correct? Who can help me solve this problem?

Problem description:

Unsteady Flow at Re = 60:Compute the unsteady solution of the scalar vorticity for a two-dimensional flow around an infinite cylinder at Re = 60. Plot the scalar vorticity after a short time interval.

The output graph is incorrect, and the value of Omega is also incorrect.

This problem has been bothering me for 2 weeks. Can anyone tell me where the problem lies? It would be best if you could help me modify the code. thank you very much indeed

function omega = flow_around_cylinder_unsteady
Re=60;
%%%%% define the grid %%%%%
n=101; m=202; % number of grid points
N=n-1; M=m-2; % number of grid intervals: 2 ghost points, theta=-h,2*pi
h=2*pi/M; % grid spacing based on theta variable
xi=(0:N)*h; theta=(-1:M)*h; % xi and theta variables on the grid
%%%%% 时间步进参数 %%%%%
t_start=0; t_end=0.5; %vortex street starts at around t=1000
tspan=[t_start t_end]; 
%%%%% 初始化涡度场 %%%%%
omega=zeros(n,m);
%%%%% 构造psi方程的矩阵A%%%%%

boundary_index = [1:n, 1:n:1+(m-1)*n, 1+(m-1)*n:m*n, n:n:n*m];
diagonals = [4 * ones(n * m, 1), -ones(n * m, 4)];
A = spdiags(diagonals, [0 -1 1 -n n], n * m, n * m);
I = speye(m * n);
A(boundary_index, :) = I(boundary_index, :);
%%%%% 查找LU分解 %%%%%
[L,U]=lu(A); clear A;
%%%%% 计算任何与时间无关的常数%%%%%


%%%%% 使用ode23进一步解决 %%%%%
options=odeset('RelTol', 1.e-03);
omega=omega(2:n-1,2:m-1); % strip boundary values for ode23
omega=omega(:); % make a column vector
[t, omega] = ode23(@(t, omega) omega_eq(omega, L, U, theta, xi, h, Re, n, m), tspan, omega, options);
%%%%% expand omega to include boundaries %%%%%
temp=zeros(n,m);
temp(2:n-1,2:m-1)=reshape(omega(end,:),n-2,m-2);
omega=temp; clear temp;
%%%%% compute stream function (needed for omega boundary values) %%%%%
omega_tilde = zeros(n, m);
omega_tilde(1, :) = 0;
omega_tilde(n, :) = exp(xi(n)) * sin(theta);
omega_tilde(:, 1) = 0;
omega_tilde(:, m) = 0;
for i = 2:n-1
    for j = 2:m-1
        omega_tilde(i, j) = h^2 * exp(2 * xi(i)) * omega(i, j);
    end
end
omega_tilde = omega_tilde(:);
psi = reshape(U \ (L \ omega_tilde), n, m);

%%%%% set omega boundary conditions %%%%%
omega(1, :) = (psi(3, :) - 8 * psi(2, :)) * 1 / (2 * h^2);
omega(n, :) = 0;
omega(:, 1) = omega(:, m-1);
omega(:, m) = omega(:, 2);

%%%%% plot scalar vorticity field %%%%%
plot_Re60(omega,t_end);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function d_omega_dt = omega_eq(omega, L, U, theta, xi, h, Re, n, m)
%%%%% expand omega to include boundary points %%%%%
temp=zeros(n,m);
index1=2:n-1; index2=2:m-1;
temp(index1,index2)=reshape(omega,n-2,m-2);
omega=temp; clear temp;
% Compute stream function
omega_tilde = zeros(n, m);
omega_tilde(1, :) = 0;
omega_tilde(n, :) = exp(xi(n)) * sin(theta);
omega_tilde(:, 1) = 0;
omega_tilde(:, m) = 0;
for i = 2:n-1
    for j = 2:m-1
        omega_tilde(i, j) = h^2 * exp(2 * xi(i)) * omega(i, j);
    end
end
omega_tilde = omega_tilde(:);
psi = reshape(U \ (L \ omega_tilde), n, m);
%%%%% compute derivatives of omega %%%%%
omega(1, :) = (psi(3, :) - 8 * psi(2, :)) * 1 / (2 * h^2);
omega(n, :) = 0;
omega(:, 1) = omega(:, m-1);
omega(:, m) = omega(:, 2);
d_omega_dt = zeros(n-2, m-2);
for i = 2:n-1
    for j = 2:m-1
        d_omega_dt(i-1, j-1) = ...
            2 * exp(-2 * xi(i)) * 1 / (h^2 * Re) * ...
            (omega(i+1, j) + omega(i-1, j) + omega(i, j+1) + omega(i, j-1) - 4 * omega(i, j)) ...
            + exp(-2 * xi(i)) / (4 * h^2) * ...
            ((psi(i+1, j) - psi(i-1, j)) * (omega(i, j+1) - omega(i, j-1)) - ...
            (psi(i, j+1) - psi(i, j-1)) * (omega(i+1, j) - omega(i-1, j)));
    end
end
d_omega_dt = d_omega_dt(:);

end

Among them, the function plot_Re60 is as follows and does not need to be modified.

function plot_Re60(omega,t_plot)
% Plot vorticity for Re = 60 from flow_past_circle_unsteady
Re=60;
% xi-theta grid
n=size(omega,1); m=size(omega,2);
N=n-1; M=m-2;
h=2*pi/M; %grid spacing
xi=(0:N)*h; theta=(-1:M)*h; %xi and theta variables on the grid
[XI, THETA] = meshgrid(xi,theta);
% x-y grid
nx=640; ny=480; %number of pixels in x and y
xmin=-1.5; xmax=10; ymax=((xmax-xmin)/2)*ny/nx; ymin=-ymax;
x=linspace(xmin,xmax,nx+1); y=linspace(ymin,ymax,ny+1);
[X,Y]=meshgrid(x,y);
%construct interpolation points
xi_i=0.5*log(X.^2+Y.^2);
theta_i=wrapTo2Pi(atan2(Y,X));
omega_xy=interp2(XI,THETA,omega',xi_i,theta_i);
%inside circle
omega_xy(xi_i<0)=0;
%set colormap for contours
levels=linspace(-1,1,1000);
v=[levels(1) levels(end)];
cmap=jet(length(levels));
cmap=flipud(cmap);
colormap(cmap);
%color contours
imagesc(x,y,omega_xy,v); hold on;
%draw black circle for cylinder
t=linspace(0,2*pi, 1000);
a=cos(t); b=sin(t);
fill(a,b,[0 0 0]);
%neaten plot
set(gca,'YDir','normal');
axis([xmin xmax ymin ymax]); set(gcf,'color','w'); axis equal; axis off;
text(xmin+0.75*(xmax-xmin),ymin+0.08*(ymax-ymin),...
['t = ', num2str(t_plot,'%3.1f')],'FontSize',22,'Color','k');
end