Combining a system of 6 equations into one variable(I) equation

Trying to combine all these six equations into one equation with variable 'I' all day, but stack with matlab error. Wanted to continue with looking at the lhs of the equation and later try to factor the coefficients of the resulting quadratic.

Here is the code:
clc
clear
% writing all six equations in terms of one equation with one variable I
syms  S E V I H R N p Lambda beta v d epsilon R0  h  r epsilon  mu sigma sigma3 s  H  eta w S0 V0

 eqn5=[h*I - epsilon*H - sigma3*H - mu*H ==0 ] ;
   H= solve(eqn5, H,'IgnoreAnalyticConstraints',true)   

 eqn6=[ r*I + epsilon*H - w*R - mu*R==0];
    R =simplifyFraction( solve(eqn6,R,'IgnoreAnalyticConstraints',true))

   eqn1=[(1-p)*Lambda-beta*I*S/N - v*S + w*R - mu*S ==0]
  S=simplifyFraction(solve(eqn1, S,'IgnoreAnalyticConstraints',true))

   eqn3=[v*S + p*Lambda - d*beta*I*V/N - mu*V==0];
    V= simplifyFraction(solve(eqn3, V,'IgnoreAnalyticConstraints',true) )

     eqn4=[eta*E - (Lambda*beta*eta*(mu +d*v-mu*p+d*mu*p))/(N*mu*(eta+mu)*(mu+v)*R0)*I ==0]
   E = simplifyFraction(solve(eqn4, E,'IgnoreAnalyticConstraints',true)) 
eqn_subs1=simplifyFraction([beta*I*S/N + d*beta*I*V/N -(eta+mu)*E == 0])
 %simlifiying the final expression for a quadratic polynomial
 eqn_subs2=expand((eqn_subs1*(mu +w)*(N*mu+I*beta*d)*(epsilon+mu+sigma3)*(I*beta+N*mu+N*v)*(mu +v)*w*mu*R0))
 % getting the LHS of the final eqution to extraction coeffcients of the
 % resuling quadratic polynomial
   LHS_eqn_subs2=lhs(eqn_subs2)