I am trying to plot differential equations using ode45 and not getting a graph

IC = [93000;1000;1000;50;50;0;500];
tspan = [0 70];
[t,Y] = ode45(@odefun, tspan, IC);
for i = 1:size(Y, 2)
    
    plot(t, Y(:, 1))
    hold on
    plot(t, Y(:, 2))
    plot(t, Y(:, 3))
    plot(t, Y(:, 4))
    plot(t, Y(:, 5))
    plot(t, Y(:, 6))
    plot(t, Y(:, 7))
    hold off
    
end
function DYdt = odefun(t,Y)
       
    S = Y(1);
    EA = Y(2);
    EB = Y(3);
    IA = Y(4);
    IS = Y(5);
    R = Y(6);
    P = Y(7);
    
    b=0.00018;
    mu=0.0004563;
    muP=0.1724;
    alpha1=0.10;
    alpha2=0.10;
    beta1=0.00414;
    beta2=0.0115;
    delta=0.7;
    O=0.7;
    psi2=0.0051;
    psi1=0.002;
    omega1=0.22;
    omega2=0.09;
    sigma=0.0018;
    gammaS=0.05;
    gammaA=0.0714;
    etaS=0.1;
    etaA=0.05;
    
    dSdt = b-((2*beta1*S*P)/(1+alpha1*P))-(2*beta2*S*(IA+IS)/1+alpha2*(IA+IS))+(psi1*EA)+(psi2*EB)-(mu*S);
    dEAdt = ((beta1*S*P)/(1+alpha1*P))+(beta2*S*(IA+IS)/1+alpha2*(IA+IS))-(psi1*EA)-(mu*EA)-(omega1*EA);
    dEBdt = ((beta1*S*P)/(1+alpha1*P))+(beta2*S*(IA+IS)/1+alpha2*(IA+IS))-(psi2*EB)-(mu*EB)-(omega2*EB);
    dIAdt = ((1-delta)*omega1*EA)-((1-O)*omega2*EB)-(mu+sigma)*IA-(gammaS*IA);
    dISdt = (delta*omega1*EA)-(O*omega2*EB)-(mu+sigma)*IS-(gammaS*IS);
    dRdt = (gammaS*IS)+(gammaA*IA)-(mu*R);
    dPdt = (etaA*IA)+(etaS*IS)-(muP*P);
    DYdt = [dSdt; dEAdt;dEBdt;dIAdt; dISdt; dRdt; dPdt];
end

Answer 1

[Torsten]

function main

IC = [93000;1000;1000;50;50;0;500];
tspan = [0 25];
[t,Y] = ode45(@odefun, tspan, IC);
plot(t, Y)   

end

function DYdt = odefun(t,Y)

S = Y(1);
EA = Y(2);
EB = Y(3);
IA = Y(4);
IS = Y(5);
R = Y(6);
P = Y(7);

b=0.00018;
mu=0.0004563;
muP=0.1724;
alpha1=0.10;
alpha2=0.10;
beta1=0.00414;
beta2=0.0115;
delta=0.7;
O=0.7;
psi2=0.0051;
psi1=0.002;
omega1=0.22;
omega2=0.09;
sigma=0.0018;
gammaS=0.05;
gammaA=0.0714;
etaS=0.1;
etaA=0.05;

dSdt = b-((2*beta1*S*P)/(1+alpha1*P))-(2*beta2*S*(IA+IS)/1+alpha2*(IA+IS))+(psi1*EA)+(psi2*EB)-(mu*S);
dEAdt = ((beta1*S*P)/(1+alpha1*P))+(beta2*S*(IA+IS)/1+alpha2*(IA+IS))-(psi1*EA)-(mu*EA)-(omega1*EA);
dEBdt = ((beta1*S*P)/(1+alpha1*P))+(beta2*S*(IA+IS)/1+alpha2*(IA+IS))-(psi2*EB)-(mu*EB)-(omega2*EB);
dIAdt = ((1-delta)*omega1*EA)-((1-O)*omega2*EB)-(mu+sigma)*IA-(gammaS*IA);
dISdt = (delta*omega1*EA)-(O*omega2*EB)-(mu+sigma)*IS-(gammaS*IS);
dRdt = (gammaS*IS)+(gammaA*IA)-(mu*R);
dPdt = (etaA*IA)+(etaS*IS)-(muP*P);
DYdt = [dSdt; dEAdt;dEBdt;dIAdt; dISdt; dRdt; dPdt];

end

[Nandita Kamath]

Thank you so much, was stuck on this one for a while!