How to plot multiple plots

My question I want to plot for different values of theta
lets 

theta= 0.1095
theta= 0.1075
theta= 0.1075
theta= 0.1055

The plot should be from 0 to 120 and the graph should start from 100 to 120

clc;clear all;close all;
% theta=0.10895;
% theta= 0.10941;
% theta=0.109428
theta= 0.1095;
YF=0.0667;
alpha= 0.29;
beta= 0.68;
gamma1=450;
gamma2=11.25;

X0=0; Y0=0.0667; Z0=0;


f=@(t,y)[-y(1)/theta+(1+alpha)*gamma1*(1-y(1))*y(2)^2+beta*gamma1*(1-y(1))*y(3)^2;...
        (YF-y(2))/theta+(1-alpha)*gamma1*(1-y(1))*y(2)^2-gamma2*y(2);...
        -y(3)/theta+beta*gamma1*(1-y(1))*y(3)^2+2*alpha*gamma1*(1-y(1))*y(2)^2-gamma2*y(3)/beta];
figure(1)
y=[X0,Y0,Z0];
[T,Y]=ode45(f,[0 5],y);
plot(T,Y(:,1),'-r','LineWidth',1.5);grid on;hold on;
plot(T,Y(:,2),'-b','LineWidth',1.5);
xlabel('Time (s)');title('\theta=0.10895')
legend('X','Y','Location','NorthEast')
figure(2)
[T,Y]=ode45(f,[100 120],y);
plot(T,Y(:,1),'-r','LineWidth',1.5);grid on;hold on;
plot(T,Y(:,2),'-b','LineWidth',1.5);
xlabel('Time (s)');title('\theta=0.10895')
legend('X','Y','Location','NorthEast')
figure(3)
[T,Y]=ode45(f,[100 120],y);
plot(Y(:,1),Y(:,2),'-b','LineWidth',1.5);grid on;hold on;
xlabel('X');ylabel('Y');title('X vs Y for \theta=0.10895')

[Rik]

Essentially you have this function:

function [T,Y1,Y2]=calculate_T_Y_Y(theta)
YF=0.0667;
alpha= 0.29;
beta= 0.68;
gamma1=450;
gamma2=11.25;
X0=0; Y0=0.0667; Z0=0;
f=@(t,y)[-y(1)/theta+(1+alpha)*gamma1*(1-y(1))*y(2)^2+beta*gamma1*(1-y(1))*y(3)^2;...
        (YF-y(2))/theta+(1-alpha)*gamma1*(1-y(1))*y(2)^2-gamma2*y(2);...
        -y(3)/theta+beta*gamma1*(1-y(1))*y(3)^2+2*alpha*gamma1*(1-y(1))*y(2)^2-gamma2*y(3)/beta];

y=[X0,Y0,Z0];
[T,Y]=ode45(f,[0 5],y);
Y1=Y(:,1);
Y2=Y(:,2);
end

So now you can calll this in a loop.

Answer 1

[Chunru]

theta_all = [0.1095,  0.1075, 0.1075, 0.1055];
for theta=theta_all
    YF=0.0667;
    alpha= 0.29;
    beta= 0.68;
    gamma1=450;
    gamma2=11.25;
    X0=0; Y0=0.0667; Z0=0;
    f=@(t,y)[-y(1)/theta+(1+alpha)*gamma1*(1-y(1))*y(2)^2+beta*gamma1*(1-y(1))*y(3)^2;...
            (YF-y(2))/theta+(1-alpha)*gamma1*(1-y(1))*y(2)^2-gamma2*y(2);...
            -y(3)/theta+beta*gamma1*(1-y(1))*y(3)^2+2*alpha*gamma1*(1-y(1))*y(2)^2-gamma2*y(3)/beta];
    figure
    y=[X0,Y0,Z0];
    [T,Y]=ode45(f,[0 5],y);
    plot(T,Y(:,1),'-r','LineWidth',1.5);grid on;hold on;
    plot(T,Y(:,2),'-b','LineWidth',1.5);
    xlabel('Time (s)');
    title(sprintf('\\theta=%f', theta))
    legend('X','Y','Location','NorthEast')
    
    figure
    [T,Y]=ode45(f,[100 120],y);
    plot(T,Y(:,1),'-r','LineWidth',1.5);grid on;hold on;
    plot(T,Y(:,2),'-b','LineWidth',1.5);
    xlabel('Time (s)');
    title(sprintf('\\theta=%f', theta))
    legend('X','Y','Location','NorthEast')

    figure
    [T,Y]=ode45(f,[100 120],y);
    plot(Y(:,1),Y(:,2),'-b','LineWidth',1.5);grid on;hold on;
    xlabel('X');ylabel('Y');
    title(sprintf('X vs Y for \\theta=%f', theta))

end