if elseif else statement not working

I have this code

r0=0.05;
k1=0.5;
k2=0.5;
K=1;
mu=0.5;
rho=0.5;
alpha=0.1;
b=0.05;
epsilon=0.25;
sigma=0.05;
eta=0.05;
pi=0.05;
omega0=0.001;

syms M m Msol positive
P(M) = epsilon*M/(1+M);
N(M) = mu+(rho*M/(1+M));
lambda(M)=1/(1+N);
r(M) = r0*(1+k1*N*(1-k2*N));
F = (m/b)*r(M)*(eta+P(M))*(1+N(M))*(1-(((alpha*sigma*m*(eta+P(M))*(1+N(M))+b*(m+alpha)*(pi*M-omega0)))/(b*alpha*sigma*K)))-(m/sigma)*(pi*M-omega0)==0;
[num,den]=numden(lhs(F));
m_num = linspace(0.01,0.9,50);
for u = 1:numel(m_num)
    Msol(u) = vpa(solve(subs(num,m,m_num(u))));
end

R0 = (b.*K.*lambda(omega0./pi))./(m_num.*(eta+P(omega0./pi)));

if R0 < 1
    l = 0;
    y = 0;
    v = 0;
else
    l = ((b.*K.*lambda(omega0./pi).*(pi.*Msol-omega0))./(alpha.*sigma.*R0.*(eta+P(omega0./pi))));
    y = ((pi.*Msol-((pi+omega0)./epsilon).*(((b.*K.*lambda(omega0./pi))./(m_num.*R0))-eta))./(sigma));
    v = ((m_num.*R0.*(eta+P(omega0./pi)).*(pi.*Msol-omega0))./(sigma.*K.*lambda(omega0./pi).*(eta+P(Msol))));
end

subplot(3,1,1)
p1 = plot (R0,l);
ylabel('l^*')
xlim([0 5]);
ylim([-0.15 0.15]);
grid on;

subplot(3,1,2)
p2 = plot (R0,y);
ylabel('y^*')
xlim([0 5]);
ylim([-0.15 0.15]);
grid on;

subplot(3,1,3)
p3 = plot (R0,v);
xlabel('R_0'), ylabel('v^*')
xlim([0 5]);
ylim([-0.15 0.15]);
grid on;

set(p1,{'LineWidth'},{1.5});
set(p2,{'LineWidth'},{1.5});
set(p3,{'LineWidth'},{1.5});

I am trying to use if,elseif,else syntax to have l,y and v zeroes for R0<1 and positive elsewhere. However, it looks like MATLAB does not recognize this part. Any help would be appreciated!

Answer 1

[Star Strider]

Use a logical vector to select the elements you want to plot, then use that as a subscript to both arguments in each plot call —

r0=0.05;
k1=0.5;
k2=0.5;
K=1;
mu=0.5;
rho=0.5;
alpha=0.1;
b=0.05;
epsilon=0.25;
sigma=0.05;
eta=0.05;
pi=0.05;
omega0=0.001;

syms M m Msol positive
P(M) = epsilon*M/(1+M);
N(M) = mu+(rho*M/(1+M));
lambda(M)=1/(1+N);
r(M) = r0*(1+k1*N*(1-k2*N));
F = (m/b)*r(M)*(eta+P(M))*(1+N(M))*(1-(((alpha*sigma*m*(eta+P(M))*(1+N(M))+b*(m+alpha)*(pi*M-omega0)))/(b*alpha*sigma*K)))-(m/sigma)*(pi*M-omega0)==0;
[num,den]=numden(lhs(F));
m_num = linspace(0.01,0.9,50);
for u = 1:numel(m_num)
    Msol(u) = vpa(solve(subs(num,m,m_num(u))));
end

R0 = (b.*K.*lambda(omega0./pi))./(m_num.*(eta+P(omega0./pi)));


l = ((b.*K.*lambda(omega0./pi).*(pi.*Msol-omega0))./(alpha.*sigma.*R0.*(eta+P(omega0./pi))));
y = ((pi.*Msol-((pi+omega0)./epsilon).*(((b.*K.*lambda(omega0./pi))./(m_num.*R0))-eta))./(sigma));
v = ((m_num.*R0.*(eta+P(omega0./pi)).*(pi.*Msol-omega0))./(sigma.*K.*lambda(omega0./pi).*(eta+P(Msol))));

Lv = R0 >= 1;               % Logical Vector

subplot(3,1,1)
p1 = plot (R0(Lv),l(Lv));
hold on
plot([0 1.02], [0 0], '-', 'Color',p1.Color, 'LineWidth',1)
hold off
ylabel('l^*')
xlim([0 5]);
ylim([-0.15 0.15]);
grid on;

subplot(3,1,2)
p2 = plot (R0(Lv),y(Lv));
hold on
plot([0 1.02], [0 0], '-', 'Color',p2.Color, 'LineWidth',1)
hold off
ylabel('y^*')
xlim([0 5]);
ylim([-0.15 0.15]);
grid on;

subplot(3,1,3)
p3 = plot (R0(Lv),v(Lv));
hold on
plot([0 1.02], [0 0], '-', 'Color',p3.Color, 'LineWidth',1)
hold off
xlabel('R_0'), ylabel('v^*')
xlim([0 5]);
ylim([-0.15 0.15]);
grid on;

set(p1,{'LineWidth'},{1.5});
set(p2,{'LineWidth'},{1.5});
set(p3,{'LineWidth'},{1.5});

EDIT — (25 Apr 2023 at 16:54)

Added constant lines from [0 1.02] in every plot.

.

[Yakoob]

Thank you Star for your help. Would it be possible for the lines l=0,y=0,v=0 to appear in the figures when R0<1? Thanks!

[Star Strider]

My pleasure!

I edited my previous code to do that.

[Yakoob]

Appreciated!!

[Star Strider]

As always, my pleasure!