Can anyone please help me with this ?

im trying to compute the unsteady solutions for the stream function and scalar vorticity for a 2d flow around an infinite cylinder at RE = 60.

This is the code i've been using:

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
%%%%% time-stepping parameters %%%%%
t_start=0; t_end=0.5; %vortex street starts at around t=1000
tspan=[t_start t_end]; 
%%%%% Initialize vorticity field %%%%%
omega=zeros(n,m);
%%%%% Construct the matrix A for psi equation %%%%%
boundary_index_bottom = 1:n;
boundary_index_left = 1:n:1+(m-1)*n;
boundary_index_top = 1+(m-1)*n:m*n;
boundary_index_right = n:n:m*n;
diagonals1 = [4*ones(n*m,1), -ones(n*m,4)];
A=spdiags(diagonals1,[0 -1 1 -n n], n*m, n*m);
I=speye(m*n);
A(boundary_index_left,:)=I(boundary_index_left,:);
A(boundary_index_right,:)=I(boundary_index_right,:);
diagonals2 = [ones(n*m,1), -ones(n*m,2)];
Aux = spdiags(diagonals2,[0 n -n],n*m,n*m);
A(boundary_index_top,:)=Aux(boundary_index_top,:);
A(boundary_index_bottom,:)=Aux(boundary_index_bottom,:);

%%%%% Find the LU decomposition %%%%%
[L,U]=lu(A); clear A;
%%%%% Compute any time-independent constants %%%%%


%%%%% advance solution using 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);
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

The code to call the function is:

omega = flow_around_cylinder_unsteady;

is there any way to fix it or improve upon it ?

[Rik]

This question looks very similar to your other question. If it is not a duplicate, feel free to reopen it and explain what is different.

[Mohamed]

This question is about unsteady solution, whereas the other question is steady solution

[per isakson]

What's the problem? What needs fixing?

I replaced plot_Re60(omega); by imagesc(omega); since I miss plot_Re60. Now the function runs and returns a result - no warnings or errors.

See TUTORIAL: how to ask a question (on Answers) and get a fast answer

[Mohamed]

the code in the gray areas cannot be changed.

i keep getting this:

this is line 132:

is there any other way to fix this ?

[Cris LaPierre]

This appears to be one of the final assignments in Mathematics for Engineers: The Capstone Course on Coursera

Answer 1

[Walter Roberson]

You have not shown enough code to be certain, but I suspect that your function plot_Re60() is defined as taking two inputs, the first of which is the array of omega values, and the second is expected to be a vector of time values. And that the problem is that at line 55 of flow_around_cylinder_unsteady, that you are not passing the time vector into plot_Re60

[Mohamed]

the code is on top.

the gray areas cannot be changed:

is it possible for you to show how to correct the code ?

i would be very thankful.

[per isakson]

Please upload the file plot_Re60.m with the paperclip icon

[Mohamed]

there is no file to upload

this is what im working with:

[Walter Roberson]

We need to see the code for the function plot_Re60() . Note the code that calls the function (you already posted that) but rather the code that defines plot_Re60

[Mohamed]

do you mean this ?

because this is all there is

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
%%%%% time-stepping parameters %%%%%
t_start=0; t_end=0.5; %vortex street starts at around t=1000
tspan=[t_start t_end]; 
%%%%% Initialize vorticity field %%%%%
omega=zeros(n,m);
%%%%% Construct the matrix A for psi equation %%%%%
boundary_index_bottom = 1:n;
boundary_index_left = 1:n:1+(m-1)*n;
boundary_index_top = 1+(m-1)*n:m*n;
boundary_index_right = n:n:m*n;
diagonals1 = [4*ones(n*m,1), -ones(n*m,4)];
A=spdiags(diagonals1,[0 -1 1 -n n], n*m, n*m);
I=speye(m*n);
A(boundary_index_left,:)=I(boundary_index_left,:);
A(boundary_index_right,:)=I(boundary_index_right,:);
diagonals2 = [ones(n*m,1), -ones(n*m,2)];
Aux = spdiags(diagonals2,[0 n -n],n*m,n*m);
A(boundary_index_top,:)=Aux(boundary_index_top,:);
A(boundary_index_bottom,:)=Aux(boundary_index_bottom,:);
%%%%% Find the LU decomposition %%%%%
[L,U]=lu(A); clear A;
%%%%% Compute any time-independent constants %%%%%
%%%%% advance solution using 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);
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

[per isakson]

Run

which plot_Re60

and show us the result

[Mohamed]

like this ?

[Walter Roberson]

We need the code of the entire function that includes the line that the error is occurring on, line 132.

[per isakson]

@Walter Roberson Sorry, I edited your comment thinking it was mine and next edited it back.

[per isakson]

@Mohamed

Yes, "like this", but I expected a different result than .

• The function, which, doesn't find plot_Re60.
• Your function, flow_around_cylinder_unsteady(), finds plot_Re60 and runs it. From the error message, , we learned that plot_Re60() encounters an error and halts the execution.

Thus, plot_Re60(), is a local function in the file, flow_around_cylinder_unsteady.m - I assume.

Please upload your file flow_around_cylinder_unsteady.m with the paperclip icon

[Mohamed]

it isnt a file

im doing this on the Mathworks Learning tool

this is all i have to work with, including the grayed out areas of the code above

[per isakson]

I'm not familiar with "Mathworks Learning tool" and I missed that you work in that environment. Furthermore, I failed to open the page, https://learningtool.mathworks.com/launch.

Somewhere under the hood there is this function, plot_Re60(), which throws the error. The background of the call, , is grey, thus you cannot test @Walter Roberson hypothesis by adding a second input argument. Or maybe, replace

%%%%% plot scalar vorticity field %%%%%

by

plot_Re60(omega,t); %%%%% plot scalar vorticity field %%%%% 

and run your function and see what happens.

This LearningTool environment doesn't it provide a way to ask questions? See https://se.mathworks.com/matlabcentral/answers/488656-zybooks-mathworks-learning-tool-submission-waiting-to-be-recorded-indefinitely

[Mohamed]

i get this

[per isakson]

@Walter Roberson assumption in his answer seems to be correct.

The statement on lines 131 and 132 cannot produce since t(end) is a scalar. I looks like the output of a vector. Something is wrong.

Two items for you:

• Didn't plot_Re60(omega); on line 55 produce any output (/error message)?
• Scroll up from line 132 until you find function ... plot_Re60( ... . Show us that line!

[Mohamed]

It's line 55 and cannot be changed

[per isakson]

That's up to you to get the algorithm right.

[per isakson]

"It's line 55 and cannot be changed" what does "it's" refer to?

I knew that line 55 cannot be changed. That's why a proposed to change line 54. Did you change line 54?

[Mohamed]

I changed it according to your instructions and got the above message I will keep trying and will let you know if I need anymore help

[per isakson]

I think you should raise the question whether the line 55 statement, plot_Re60(omega), is correct. But I don't know with whom. You instructor?

One cannot be sure that the trick of changing line 54 doesn't have some serious side effects.

[Hugo]

I am also struggling with this assignment. The code for plot_Re60 is the following:

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

[Walter Roberson]

The call to plot_Re60 (which cannot be changed) only passes in a single parameter, but it appears from @Hugo that the function expects two parameters to be passed, with the second one named t_plot inside the function.

The code posted by Hugo uses t_plot only once, right near the bottom, in the lines

text(xmin+0.75*(xmax-xmin),ymin+0.08*(ymax-ymin),...
['t = ', num2str(t_plot,'%3.1f')],'FontSize',22,'Color','k');

If only a single parameter was being passed to the function, then Yes, you would have a problem: you would get an error message about not enough input parameters.

The code posted at https://www.mathworks.com/matlabcentral/answers/1815940-can-anyone-please-help-me-with-this#comment_2393145 by @Mohamed shows very similar code, but does not use a variable named t_plot at that point in the code, and instead uses t . If the only change between Mohamed's code and Hugo's code was that for Mohamed, the parameter was named t instead of t_plot then there would be the point that in the code posted by Hugo,

t=linspace(0,2*pi, 1000);

and so if that statement existed in Mohamed's code but the text() call referred to t then there would not be any error about insufficient parameters -- since the missing parameter would have been overwritten by the assignment of the linspace result.

We deduce that the function code for Mohamed is not the same as the function code for Hugo.

If a time is needed by the plotting function, but that it is not permitted to change the line that calls the function, but it is permitted to change the function, then I can see two possibilities at the moment:

1. You could alter the plotting function to not expect a second parameter, and instead to use evalin('caller') to read out the time from the calling function. This is not generally a good practice, but if the teacher refuses to fix the calling code then you have to hack; OR
2. since you are permitted to modify the comment immediately above the call to the plotting function, you could insert the proper call there. And then you could code the plotting function to detect that nargin < 2 and just to return in that case. So the badly-coded call to the plotting function would have no effect on the output, and the proper call to the plotting function would do the work.

Answer 2

[Hugo]

@Mohamed did you end up managing to finish the assignment with the correct code? So yes, could you please share the code with me?

Kind Regards,

Hugo

[Basant Ale]

Hello there, did you finish the assignment? If yes, can you help me out because I am facing the same problem?

[Marco Aguilar Condori]

Excuse me, were you able to solve it? I'm really stuck on this problem and I can't figure it out

[Walter Roberson]

@Marco Aguilar Condori see my update at https://www.mathworks.com/matlabcentral/answers/1815940-can-anyone-please-help-me-with-this#comment_2905049 just posted.

Of course better would be to get the learning system to repair the bug in the code.

[Pravar Agnihotri]

Hi! Unfortunately I am still stuck, I have a code similar to Mohammed's one above. Would be awesome if you could show what changes to make. The 2nd option in your update did not work :(

[Karthik Gk]

Hi I am Also Unable to complete this Assignment . Please any one kind soul please help me to complete this problem . i would be really thankfull. please.....

#Its showing T=0.00.00 but i think it has to show like t=0.00 only pleae help me

flow_around_cylinder_unsteady();

Unrecognized function or variable 'plot_Re60'.

Error in solution>flow_around_cylinder_unsteady (line 87)
plot_Re60(omega, 0); % Plot the scalar vorticity field at t = 0

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
%%%%% time-stepping parameters %%%%%
t_start=0; t_end=0.5; %vortex street starts at around t=1000
tspan=[t_start t_end]; 
 %%%%% plot scalar vorticity field %%%%% 
omega=zeros(n,m);
%%%%% Construct the matrix A for psi equation %%%%%
boundary_index_bottom = 1:n;
boundary_index_left = 1:n:1 + (m - 1) * n;
boundary_index_top = 1 + (m - 1) * n:m * n;
boundary_index_right = n:n:m * n;
diagonals1 = [4 * ones(n * m, 1), -ones(n * m, 4)];
A = spdiags(diagonals1, [0 -1 1 -n n], n * m, n * m);
I = speye(m * n);
A(boundary_index_left, :) = I(boundary_index_left, :);
A(boundary_index_right, :) = I(boundary_index_right, :);
diagonals2 = [ones(n * m, 1), -ones(n * m, 2)];
Aux = spdiags(diagonals2, [0 n -n], n * m, n * m);
A(boundary_index_top, :) = Aux(boundary_index_top, :);
A(boundary_index_bottom, :) = Aux(boundary_index_bottom, :);

%%%%% Find the LU decomposition %%%%%
[L,U]=lu(A); clear A;
%%%%% Compute any time-independent constants %%%%%


%%%%% advance solution using 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
function d_omega_dt = omega_eq(t, omega, L, U, theta, xi, h, Re, n, m)
    % Initialize the derivative vector
    d_omega_dt = zeros(n-2, m-2);

    % Check if t is 0
    if t == 0
        % Compute derivatives of omega only at t = 0
        % Your computation code here
    else
        % For other time values, return empty derivatives
        d_omega_dt = zeros(n-2, m-2);
    end

    % Reshape d_omega_dt into a column vector
    d_omega_dt = d_omega_dt(:);
end


[t, omega] = ode23(@(t, omega) omega_eq(t, 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);

% Limit the value of t to 0.00
t = 0.00;
plot_Re60(omega, 0); % Plot the scalar vorticity field at t = 0


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

[Walter Roberson]

Inserting the plot routine but otherwise not changing the code:

flow_around_cylinder_unsteady();

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
%%%%% time-stepping parameters %%%%%
t_start=0; t_end=0.5; %vortex street starts at around t=1000
tspan=[t_start t_end]; 
 %%%%% plot scalar vorticity field %%%%% 
omega=zeros(n,m);
%%%%% Construct the matrix A for psi equation %%%%%
boundary_index_bottom = 1:n;
boundary_index_left = 1:n:1 + (m - 1) * n;
boundary_index_top = 1 + (m - 1) * n:m * n;
boundary_index_right = n:n:m * n;
diagonals1 = [4 * ones(n * m, 1), -ones(n * m, 4)];
A = spdiags(diagonals1, [0 -1 1 -n n], n * m, n * m);
I = speye(m * n);
A(boundary_index_left, :) = I(boundary_index_left, :);
A(boundary_index_right, :) = I(boundary_index_right, :);
diagonals2 = [ones(n * m, 1), -ones(n * m, 2)];
Aux = spdiags(diagonals2, [0 n -n], n * m, n * m);
A(boundary_index_top, :) = Aux(boundary_index_top, :);
A(boundary_index_bottom, :) = Aux(boundary_index_bottom, :);

%%%%% Find the LU decomposition %%%%%
[L,U]=lu(A); clear A;
%%%%% Compute any time-independent constants %%%%%


%%%%% advance solution using 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
function d_omega_dt = omega_eq(t, omega, L, U, theta, xi, h, Re, n, m)
    % Initialize the derivative vector
    d_omega_dt = zeros(n-2, m-2);

    % Check if t is 0
    if t == 0
        % Compute derivatives of omega only at t = 0
        % Your computation code here
    else
        % For other time values, return empty derivatives
        d_omega_dt = zeros(n-2, m-2);
    end

    % Reshape d_omega_dt into a column vector
    d_omega_dt = d_omega_dt(:);
end


[t, omega] = ode23(@(t, omega) omega_eq(t, 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);

% Limit the value of t to 0.00
t = 0.00;
plot_Re60(omega, 0); % Plot the scalar vorticity field at t = 0


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



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

[Walter Roberson]

I guess I do not understand what is being asked.