Trying to solve system of ODE

edit

The given code is used to find out the solution to the system of equations. All equations need to be equated to 0 and soln must be found.

However, when I run this code I am not getting any soln.

var('S,E1,E2,V,Ip,Iexp,T,R,p1,p2,p3,beta,delta,theta,mu,w1,w2,gammap, gammaexp, sigma,m,r, q, n, i, j ');
N = S+V+E1+E2+Ip+Iexp+R;
eqS = pi-(mu*S)-(p1*beta*S*(Ip+Iexp))-(p2*beta*S*(Ip+Iexp))-(p3*m*S)+(theta*V);
eqE1 = (p1*beta*S*(Ip+Iexp))-(mu*E1)-(w1*E1)+((1-n)*q*beta*V*(Ip+Iexp))+((1-j)*i*beta*R*(Ip+Iexp));
eqE2 = (p2*beta*S*(Ip+Iexp))-(mu*E2)-(w2*E2)+(n*q*beta*V*(Ip+Iexp))+(j*i*beta*R*(Ip+Iexp));
eqV=(p3*m*S)-(theta*V)-(mu*V)-((1-n)*q*beta*V*(Ip+Iexp))-(n*q*beta*V*(Ip+Iexp));
eqIp = (delta*w1*E1)+(delta*w2*E2)-(gammap*Ip)-((mu+sigma)*Ip);    
eqIexp = ((1-delta)*w1*E1)+((1-delta)*w2*E2)-((mu+sigma)*Iexp)-(gammaexp*Iexp);    
eqR=(gammap*Ip)+(gammaexp*Iexp)-(mu*R)-((1-j)*i*beta*R*(Ip+Iexp))-(j*i*beta*R*(Ip+Iexp));    
pretty_print((eqS+eqE1+eqE2+eqV+eqIp+eqIexp+eqR).full_simplify())    
sistema = [eqS==0, eqE1 == 0, eqE2==0, eqV==0,  eqIp == 0,eqIexp==0, eqR ];    
sol = solve(sistema, S,E1, E2,V, Ip,Iexp, R); pretty_print(sol)    
DFE = sol[2]; pretty_print(DFE)    
EE = sol[0]; pretty_print(EE)