"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
