I'm not sure why you would use a Newton method to solve what solve will do directly, in an apparent algebraic form.
F(1) = p*M - beta1*S*E - beta2*S*I - (rho + mu)*S;
F(2) = beta1*S*E + beta2*S*I - (omega + alpha1 + mu)*E;
F(3) = omega*E - (gamma + alpha2 + mu)*I;
F(4) = alpha1*E + alpha2*I - (theta + mu)*Q;
F(5) = (1 - p)*M + rho*S + gamma*I + theta*Q - mu*V;
sol = solve(F,[S E I Q V])
sol =
S: [2x1 sym]
E: [2x1 sym]
I: [2x1 sym]
Q: [2x1 sym]
V: [2x1 sym]
[sol.S,sol.E,sol.I,sol.Q,sol.V]
The first solution has three of the unknowns as identically zero, with only S and V non-zero. That is probably a degenerate solution that you don't care to find. But the second solution is non-trivial.
If I look at the equations carefully, It seems there are only two equations that are not purely linear in the parameter. The first two. And there we see you have the terms S*E and S*I as products.
F(1) = p*M - beta1*S*E - beta2*S*I - (rho + mu)*S;
F(2) = beta1*S*E + beta2*S*I - (omega + alpha1 + mu)*E;
But, if we just add those first two equations, so form F(1) + F(2), the products go completely away. My guess is this makes the problem now almost solvable by hand. In the end, it will reduce to a quadratic, and so, you end up with two roots.