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.
.
