Roots of an equation when two parameters are changed

I have this equation

r0=0.05;
k1=0.5; 
k2=0.5;
mu=0.5;
rho=0.5;
epsilon=0.25; 
K=1; 
alpha=0.1; 
q=0.1;
b=0.8; 
zeta=0.075;
omega0=0.001

syms M sigma eta Msol positive

Nut(M) = mu+(rho*M/(1+M));
gro(M) = r0*(1+k1*Nut(M)*(1-k2*Nut(M)));
lam(M) = 1/(1+Nut(M));
P(M) = epsilon*M/(1+M);
eqn = (q./b).*gro(M).*(eta+P(M)).*(1+Nut(M)).*(1-(((alpha.*sigma.*q.*(eta+P(M)).*(1+Nut(M))+b.*(q+alpha).*(zeta.*M-omega0)))./(b.*alpha.*sigma.*K)))-(q./sigma).*(zeta.*M-omega0)==0;

I can compute the roots of this equation when sigma is varied and eta = 0.05 using the command

[num,den]=numden(lhs(eqn));
sigma_num = linspace(0.01,1,20);
for u = 1:numel(sigma_num)
    Msol(u) = vpa(solve(subs(num,sigma,sigma_num(u))));
end

which will produce 20 roots, each one is related to each value of sigma. Now, I would like to compute the roots of the equation when two parameters are varied, say sigma and eta such that eta_num = linspace(0.01,1,20) as well. I believe I should end up having 400 roots but How I can do that? Any help is appreciated! Many thanks!

Answer 1

[Torsten]

"num" is a polynomial of degree 7 in M. Thus it has seven roots. But your code only gives one of these seven. So I hope you get the single root out of seven that you are after.

r0=0.05;
k1=0.5; 
k2=0.5;
mu=0.5;
rho=0.5;
epsilon=0.25; 
K=1; 
alpha=0.1; 
q=0.1;
b=0.8; 
zeta=0.075;
omega0=0.001;
syms M sigma eta Msol positive
Nut(M) = mu+(rho*M/(1+M));
gro(M) = r0*(1+k1*Nut(M)*(1-k2*Nut(M)));
lam(M) = 1/(1+Nut(M));
P(M) = epsilon*M/(1+M);
eqn = (q./b).*gro(M).*(eta+P(M)).*(1+Nut(M)).*(1-(((alpha.*sigma.*q.*(eta+P(M)).*(1+Nut(M))+b.*(q+alpha).*(zeta.*M-omega0)))./(b.*alpha.*sigma.*K)))-(q./sigma).*(zeta.*M-omega0)==0;
[num,den]=numden(lhs(eqn));
sigma_num = linspace(0.01,1,20);
eta_num = linspace(0.01,1,20);
n = numel(eta_num);
m = numel(sigma_num);
M_num = zeros(n,m);
%M_num = cell(n,m);
for i = 1:n
    for j = 1:m
        M_num(i,j) = vpa(solve(subs(num,[eta sigma],[eta_num(i),sigma_num(j)])==0,M));
        %M_num{i,j} = roots(sym2poly(subs(num,[eta sigma],[eta_num(i),sigma_num(j)])));
    end
end
M_num

M_num = 20×20

0.0135 0.0146 0.0157 0.0168 0.0180 0.0193 0.0206 0.0219 0.0233 0.0248 0.0263 0.0279 0.0296 0.0314 0.0333 0.0352 0.0373 0.0395 0.0418 0.0442 0.0143 0.0194 0.0247 0.0302 0.0360 0.0420 0.0482 0.0547 0.0615 0.0686 0.0760 0.0837 0.0918 0.1002 0.1090 0.1181 0.1276 0.1376 0.1479 0.1587 0.0150 0.0240 0.0333 0.0431 0.0532 0.0638 0.0748 0.0862 0.0982 0.1106 0.1235 0.1369 0.1509 0.1654 0.1804 0.1959 0.2120 0.2286 0.2458 0.2635 0.0157 0.0284 0.0416 0.0554 0.0697 0.0847 0.1003 0.1164 0.1332 0.1507 0.1687 0.1874 0.2068 0.2268 0.2474 0.2687 0.2905 0.3130 0.3360 0.3597 0.0164 0.0326 0.0495 0.0672 0.0856 0.1047 0.1246 0.1452 0.1666 0.1887 0.2116 0.2352 0.2595 0.2845 0.3102 0.3366 0.3636 0.3912 0.4194 0.4481 0.0170 0.0366 0.0571 0.0785 0.1007 0.1238 0.1478 0.1727 0.1983 0.2249 0.2522 0.2803 0.3091 0.3387 0.3690 0.4000 0.4315 0.4637 0.4965 0.5297 0.0176 0.0405 0.0644 0.0893 0.1152 0.1421 0.1699 0.1987 0.2284 0.2590 0.2905 0.3227 0.3558 0.3895 0.4240 0.4591 0.4948 0.5311 0.5679 0.6053 0.0182 0.0442 0.0713 0.0996 0.1289 0.1594 0.1909 0.2234 0.2569 0.2913 0.3266 0.3627 0.3995 0.4371 0.4754 0.5142 0.5537 0.5937 0.6342 0.6752 0.0187 0.0477 0.0779 0.1094 0.1420 0.1759 0.2108 0.2468 0.2838 0.3217 0.3605 0.4001 0.4405 0.4816 0.5233 0.5656 0.6085 0.6519 0.6958 0.7401 0.0193 0.0510 0.0842 0.1187 0.1545 0.1915 0.2296 0.2689 0.3091 0.3503 0.3924 0.4353 0.4789 0.5232 0.5681 0.6135 0.6595 0.7060 0.7529 0.8002

[Fares]

Thank you Torsten!