How do I write a for loop for different initial values of omega and then plot multiple lines?

I'm running a program to simulate a simple pendulum. I am plotting a phase space diagram but I would like to have multiple different initial values of omega plotted, and my professor's advice was to use a for loop instead of plugging in many different values. Can anyone show me how to do this? I have very few weeks of experience using Matlab.

%Simple Pendulum using Euler Method
clear all; help pendul;

%Initial Values
theta0 = input('Enter initial angle (in degrees): ');
theta = theta0*pi/180; %Convert angle to radians
omega = 0;

%Physical Constants
g_over_L = 1;   %The constant g/L
time = 0;       %Initial time
irev = 0;       %Used to count the number of reversals
tau = input('Enter time step: ');

%Loop over desired number of steps with given time step
nstep = input('Enter number of time steps: ');
for istep=1:nstep
    
    %Record omega and theta for phase space plotting
    displace_plot(istep) = theta;
    omega_plot(istep) = omega;
    
    %Compute new position and velocity
    accel = -g_over_L*sin(theta); %Gravitational Acceleration
    theta_old = theta;          %Save Previous Angle
    theta = theta + tau*omega; 
    omega = omega + tau*accel;
end

%Graph omega versus theta
figure(2)
plot(displace_plot,omega_plot, '+');
xlim([-pi,pi]);
xlabel('\theta (rads)'); ylabel('\omega (rad/s)');

[VBBV]

Code with multiple omega needs a vector along with variable inside for loop

 % if true
    % code
 % end
%Simple Pendulum using Euler Method
clear all; help pendul;
%Initial Values
theta0 = input('Enter initial angle (in degrees): ');
theta = theta0*pi/180; %Convert angle to radians
omega = 0:1:20; % speed in rpm
omega = 2*pi.*omega/60;
%Physical Constants
g_over_L = 1;   %The constant g/L
time = 0;       %Initial time
irev = 0;       %Used to count the number of reversals
tau = input('Enter time step: ');
%Loop over desired number of steps with given time step
nstep = input('Enter number of time steps: ');
for j = 1:length(omega)
Om = omega(j); % assign each omega
for istep=1:nstep
    %Record omega and theta for phase space plotting
    displace_plot(istep,j) = theta;
    omega_plot(istep,j) = Om;

    %Compute new position and velocity
    accel = -g_over_L*sin(theta); %Gravitational Acceleration
    theta_old = theta;          %Save Previous Angle
    theta = theta + tau*Om; % increment omega with angle
    Om = Om+ tau*accel; % increment omega with accel
end
clear Om;
end
%Graph omega versus theta
figure(2)
plot(displace_plot,omega_plot);
%xlim([-pi,pi]);
xlabel('\theta (rads)'); ylabel('\omega (rad/s)');

[Nicholas Davis]

This works very well, Vasishta, and I thank you for that. What addition can I make to the code to make the omega run over all integers, much like in the image attached. It only seems to run for the positive integers and a small change makes it run for only the negative integers. I'd like it to run over both. Thanks!