(* -------------- Majority Rule Group Size G odd -------------- *)
(* Parameters setting and initial conditions *)
G = 3;
NS = G-1;
(Subscript[ρ, A])=1;
α=((Subscript[ρ, A])*g)^-1;
β=((Subscript[ρ, B])*g)^-1;

(* pXY:probability that an agent X switch opinion towards Y as a result of the Majority Rule group size G=NS+1 *)
(Subscript[p, A]):=((Subscript[d, A])/(((Subscript[d, A])+(Subscript[d, B]))));
(Subscript[p, B]):=((Subscript[d, B])/(((Subscript[d, A])+(Subscript[d, B]))));
(* pAA :=Sum[NS!/(i!*(NS-i)!)*pA^i*pB^(NS-i),{i,Floor[NS/2],NS}];*)
(Subscript[p, BA]):=(∑_(i=0)^(Floor[(NS/2)]-1))((NS!)/(i!*(NS-i)!))*((Subscript[p, A])^((NS-i)))*((Subscript[p, B])^i);
(Subscript[p, AB]):=(∑_(i=0)^(Floor[(NS/2)]-1))((NS!)/(i!*(NS-i)!))*((Subscript[p, A])^i)*((Subscript[p, B])^((NS-i)));
(Subscript[p, AA]):=1-(Subscript[p, AB]);
(Subscript[p, BB]):=1-(Subscript[p, BA]);

(* Definition of the Steady State System *)
steadySystem={
σ*(Subscript[e, A])-α*(Subscript[d, A])==0,
σ*(Subscript[e, B])-β*(Subscript[d, B])==0,
(Subscript[p, AA])*α*(Subscript[d, A])+(Subscript[p, BA])*β*(Subscript[d, B])-σ*(Subscript[e, A])==0,
(Subscript[p, AB])*α*(Subscript[d, A])+(Subscript[p, BB])*β*(Subscript[d, B])-σ*(Subscript[e, B])==0};
conditions={
(Subscript[d, A])+(Subscript[d, B])+(Subscript[e, A])+(Subscript[e, B])-1==0,
(Subscript[d, A])≥ 0,
(Subscript[d, B])≥ 0,
(Subscript[e, A])≥ 0,
(Subscript[e, B])≥ 0,
σ>0,
g>0,
1>(Subscript[ρ, B])>0};
variables ={(Subscript[d, A]),(Subscript[d, B]),(Subscript[e, A]),(Subscript[e, B])};

(* Compute steady states *)
rule=ConditionalExpression[a_,b__]->a;
equilibria=Simplify[Solve[Rationalize[Join[steadySystem,conditions]],variables]/.rule];
Print["There are ",Length[equilibria]," steady states"];
Print[equilibria];

(* Compute eigenvalues *)
steadySystemLHS=steadySystem/.Equal[x_,y_]->x;
jac=Outer[D,steadySystemLHS,variables];
eigenValues=Map[Eigenvalues[jac/.#]&,equilibria];

(* Stability analysis *)
For[eq=1,eq<Length[eigenValues]+1,eq++,
ev=eigenValues[[eq]];
Print["Equilibrium[",eq,"]:",equilibria[[eq]]];
Print["Eigenvalues:",ev];
For[evIt=2,evIt<Length[ev]+1,evIt++,
eval=Rationalize[ev[[evIt]]];
conditions={FullSimplify[eval]<0,σ>0,g>0,1>(Subscript[ρ, B])>0};
variables={σ, g,(Subscript[ρ, B])};
red=Reduce[conditions,variables];
Print["Condition satisfied(ev[",evIt,"]<0): ",red];];
Print["--------------------------------------"];
];