(* Proof of Lemma 1 Gibbs (2015); Last Updated 25/09/2015 *)
(* The Tmap *)
T[M1_List,M2_List,w_List]:=Module[{A,B,CC,Ma,Mb,TM1,TM2,TW1,TW2,theta, thetaj,Ez,u1,u2,rho},
A=μ+α (w[[1]]M1[[1]]+w[[2]]M2[[1]]);
B=α w[[1]]M1[[2]]ρ1+ζ1;
CC=α w[[2]] M2[[2]]ρ2+ζ2;
Ma={{1,0},{0,z1}}//Inverse;
Mb={{1,0},{0,z2}}//Inverse;

TM1=Ma.{{A},{B z1+CC s}};
TM2=Mb.{{A},{B s+CC z2}};
theta ={{M1[[1]]},{M1[[2]]}};
thetaj ={{M2[[1]]},{M2[[2]]}};
Ez ={{1,0,0},{0,z1,s},{0,s,z2}};
u1 ={{1,0,0},{0,1,0}};
u2 ={{1,0,0},{0,0,1}};
rho ={{1,0,0},{0,ρ1,0},{0,0,ρ2}};
TW1 = Inverse[-Transpose[theta].u1.rho.Ez.Transpose[rho].(α*rho-IdentityMatrix[3]).Transpose[u1].theta]*(Transpose[theta].u1.rho.Ez.rho.{{μ},{ζ1},{ζ2}} + w[[2]]*Transpose[theta].u1.rho.Ez.Transpose[rho].(α*rho-IdentityMatrix[3]).Transpose[u2].thetaj);
TW2 = Inverse[-Transpose[thetaj].u2.rho.Ez.Transpose[rho].(α*rho-IdentityMatrix[3]).Transpose[u2].thetaj]*(Transpose[thetaj].u2.rho.Ez.rho.{{μ},{ζ1},{ζ2}} + w[[1]]*Transpose[thetaj].u2.rho.Ez.Transpose[rho].(α*rho-IdentityMatrix[3]).Transpose[u1].theta);
{TM1[[1]]-M1[[1]],TM1[[2]]-M1[[2]],TM2[[1]]-M2[[1]],TM2[[2]]-M2[[2]],Flatten@TW1-w[[1]], Flatten@TW2-w[[2]]}



];
s=0;
μ=0;

(* Transforming the system into normal form. See page 35 of Gibbs (2015)

TransFormSys[a1_,a2_,b1_,b2_,γ1_,γ2_,α_]:=Module[{cc,ss1,ss2,ss3,ss4,ss5},
cc=a1 (-1+α γ1)+a2 α γ2;
ss1=a1 α γ1+a2 (-1+α γ2);
ss2=ζ1+b1 (-1+α γ1 ρ1);
ss3=ζ2+b2 (-1+α γ2 ρ2);
ss4=-γ1+((a1 a2 (-1+α) γ2+b1 z1 ζ1 ((ρ1)^2))/(-(a1^2) (-1+α)+(b1^2) z1 ((ρ1)^2) (1-α ρ1)));
ss5=-γ2+((a1 a2 (-1+α) γ1+b2 z2 ζ2 ((ρ2)^2))/(-(a2^2) (-1+α)+(b2^2) z2 ((ρ2)^2) (1-α ρ2)));
{cc,ss1,ss2,ss3,ss4,ss5}
];
TransFormSys[u1,u2,v1+ζ1/(1-(ϵ+1/2)ρ1),v2+ζ2/(1-(ϵ+1/2)ρ2
),w1+1,w2+1,ϵ+1/2]/.{u1->0,u2->0,v1->0,v2->0,w1->0,w2->0,ϵ->0}//Simplify
NewSystem[u1_,u2_,v1_,v2_,w1_,w2_,ϵ_]:=TransFormSys[u1,u2,v1+ζ1/(1-(ϵ+1/2)ρ1),v2+ζ2/(1-(ϵ+1/2)ρ2),w1+1,w2+1,ϵ+1/2];
NewSystem[u1,u2,v1,v2,w1,w2,ϵ]//Simplify //MatrixForm

DTI:=D[Flatten@NewSystem[u1,u2,v1,v2,w1,w2,ϵ],{{u1,u2,v1,v2,w1,w2}}]/.{u1->0,u2->0,v1->0,v2->0,w1->0,w2->0,ϵ->0}

V:=JordanDecomposition[DTI][[1]]
EigDTI:=JordanDecomposition[DTI][[2]]

EigDTI

Expand[NewSystem[u1,u2,v1,v2,w1,w2,ϵ][[5]]]

Order2Terms[u1_,u2_,v1_,v2_,w1_,w2_]:=Module[{f1,f2,f3,f4,f5,f6},
f1=((u1 w1)/2)+((u2 w2)/2)+u1 ϵ+u2 ϵ+u1 w1 ϵ+u2 w2 ϵ;
f2=((u1 w1)/2)+((u2 w2)/2)+u1 ϵ+u2 ϵ+u1 w1 ϵ+u2 w2 ϵ;
f3=((v1 w1 ρ1)/2)+v1 ϵ ρ1+v1 w1 ϵ ρ1-((ζ1)/(1-((1/2)+ϵ) ρ1))+((ζ1 ρ1)/(2 (1-((1/2)+ϵ) ρ1)))+((w1 ζ1 ρ1)/(2 (1-((1/2)+ϵ) ρ1)))+((ϵ ζ1 ρ1)/(1-((1/2)+ϵ) ρ1))+((w1 ϵ ζ1 ρ1)/(1-((1/2)+ϵ) ρ1));
f4=((v2 w2 ρ2)/2)+v2 ϵ ρ2+v2 w2 ϵ ρ2-((ζ2)/(1-((1/2)+ϵ) ρ2))+((ζ2 ρ2)/(2 (1-((1/2)+ϵ) ρ2)))+((w2 ζ2 ρ2)/(2 (1-((1/2)+ϵ) ρ2)))+((ϵ ζ2 ρ2)/(1-((1/2)+ϵ) ρ2))+((w2 ϵ ζ2 ρ2)/(1-((1/2)+ϵ) ρ2));
f5=-((u1 u2)/(2 (-(u1^2) (-(1/2)+ϵ)+z1 ((ρ1)^2) (1-((1/2)+ϵ) ρ1) (((v1+((ζ1)/(1-((1/2)+ϵ) ρ1))))^2))))-((u1 u2 w2)/(2 (-(u1^2) (-(1/2)+ϵ)+z1 ((ρ1)^2) (1-((1/2)+ϵ) ρ1) (((v1+((ζ1)/(1-((1/2)+ϵ) ρ1))))^2))))+((u1 u2 ϵ)/(-(u1^2) (-(1/2)+ϵ)+z1 ((ρ1)^2) (1-((1/2)+ϵ) ρ1) (((v1+((ζ1)/(1-((1/2)+ϵ) ρ1))))^2)))+((u1 u2 w2 ϵ)/(-(u1^2) (-(1/2)+ϵ)+z1 ((ρ1)^2) (1-((1/2)+ϵ) ρ1) (((v1+((ζ1)/(1-((1/2)+ϵ) ρ1))))^2)))+((v1 z1 ζ1 ((ρ1)^2))/(-(u1^2) (-(1/2)+ϵ)+z1 ((ρ1)^2) (1-((1/2)+ϵ) ρ1) (((v1+((ζ1)/(1-((1/2)+ϵ) ρ1))))^2)))+((z1 ((ζ1)^2) ((ρ1)^2))/((1-((1/2)+ϵ) ρ1) (-(u1^2) (-(1/2)+ϵ)+z1 ((ρ1)^2) (1-((1/2)+ϵ) ρ1) (((v1+((ζ1)/(1-((1/2)+ϵ) ρ1))))^2))));
f6=-((u1 u2)/(2 (-(u2^2) (-(1/2)+ϵ)+z2 ((ρ2)^2) (1-((1/2)+ϵ) ρ2) (((v2+((ζ2)/(1-((1/2)+ϵ) ρ2))))^2))))-((u1 u2 w1)/(2 (-(u2^2) (-(1/2)+ϵ)+z2 ((ρ2)^2) (1-((1/2)+ϵ) ρ2) (((v2+((ζ2)/(1-((1/2)+ϵ) ρ2))))^2))))+((u1 u2 ϵ)/(-(u2^2) (-(1/2)+ϵ)+z2 ((ρ2)^2) (1-((1/2)+ϵ) ρ2) (((v2+((ζ2)/(1-((1/2)+ϵ) ρ2))))^2)))+((u1 u2 w1 ϵ)/(-(u2^2) (-(1/2)+ϵ)+z2 ((ρ2)^2) (1-((1/2)+ϵ) ρ2) (((v2+((ζ2)/(1-((1/2)+ϵ) ρ2))))^2)))+((v2 z2 ζ2 ((ρ2)^2))/(-(u2^2) (-(1/2)+ϵ)+z2 ((ρ2)^2) (1-((1/2)+ϵ) ρ2) (((v2+((ζ2)/(1-((1/2)+ϵ) ρ2))))^2)))+((z2 ((ζ2)^2) ((ρ2)^2))/((1-((1/2)+ϵ) ρ2) (-(u2^2) (-(1/2)+ϵ)+z2 ((ρ2)^2) (1-((1/2)+ϵ) ρ2) (((v2+((ζ2)/(1-((1/2)+ϵ) ρ2))))^2))));
{f1,f2,f3,f4,f5,f6}
];

V.{{u1},{u2},{v1},{v2},{w1},{w2}}

EigDTI.{{u1},{u2},{v1},{v2},{w1},{w2}}

NSystem[u1_,u2_,v1_,v2_,w1_,w2_]:=Inverse[V].Order2Terms[-u1+u2,u1+u2,-((v1 ζ1 ((Sqrt[-8+ρ1]) (Sqrt[ρ1])-ρ1))/(2 (-2+ρ1)))+((v2 ζ1 ((Sqrt[-8+ρ1]) (Sqrt[ρ1])+ρ1))/(2 (-2+ρ1))),-((w1 ζ2 ((Sqrt[-8+ρ2]) (Sqrt[ρ2])-ρ2))/(2 (-2+ρ2)))+((w2 ζ2 ((Sqrt[-8+ρ2]) (Sqrt[ρ2])+ρ2))/(2 (-2+ρ2))),v1+v2,w1+w2]

(* Equation 37 on page 35 of Gibbs (2015) *)

N1[nw2,ϵ]:={{D[h1[nw2,ϵ],nw2]*NSystem[h1[nw2,ϵ],nw2,h2[nw2,ϵ],h3[nw2,ϵ],h4[nw2,ϵ],h5[nw2,ϵ]][[2]]-(-1)*h1[nw2,ϵ]-NSystem[h1[nw2,ϵ],nw2,h2[nw2,ϵ],h3[nw2,ϵ],h4[nw2,ϵ],h5[nw2,ϵ]][[1]]},{D[h2[nw2,ϵ],nw2]*NSystem[h1[nw2,ϵ],nw2,h2[nw2,ϵ],h3[nw2,ϵ],h4[nw2,ϵ],h5[nw2,ϵ]][[2]]-(1/4) (-4-(Sqrt[-8+ρ1]) (Sqrt[ρ1])+ρ1)*h2[nw2,ϵ]-NSystem[h1[nw2,ϵ],nw2,h2[nw2,ϵ],h3[nw2,ϵ],h4[nw2,ϵ],h5[nw2,ϵ]][[3]]},{D[h3[nw2,ϵ],nw2]*NSystem[h1[nw2,ϵ],nw2,h2[nw2,ϵ],h3[nw2,ϵ],h4[nw2,ϵ],h5[nw2,ϵ]][[2]]-(1/4) (-4+(Sqrt[-8+ρ1]) (Sqrt[ρ1])+ρ1)*h3[nw2,ϵ]-NSystem[h1[nw2,ϵ],nw2,h2[nw2,ϵ],h3[nw2,ϵ],h4[nw2,ϵ],h5[nw2,ϵ]][[4]]},{D[h4[nw2,ϵ],nw2]*NSystem[h1[nw2,ϵ],nw2,h2[nw2,ϵ],h3[nw2,ϵ],h4[nw2,ϵ],h5[nw2,ϵ]][[2]]-(1/4) (-4-(Sqrt[-8+ρ2]) (Sqrt[ρ2])+ρ2)*h4[nw2,ϵ]-NSystem[h1[nw2,ϵ],nw2,h2[nw2,ϵ],h3[nw2,ϵ],h4[nw2,ϵ],h5[nw2,ϵ]][[5]]},{D[h5[nw2,ϵ],nw2]*NSystem[h1[nw2,ϵ],nw2,h2[nw2,ϵ],h3[nw2,ϵ],h4[nw2,ϵ],h5[nw2,ϵ]][[2]]-(1/4) (-4+(Sqrt[-8+ρ2]) (Sqrt[ρ2])+ρ2)*h5[nw2,ϵ]-NSystem[h1[nw2,ϵ],nw2,h2[nw2,ϵ],h3[nw2,ϵ],h4[nw2,ϵ],h5[nw2,ϵ]][[6]]}}

(* Taylor Approximation of the center manifold *)

Approx[t_]:=Module[{rule,hh1,hh2,hh3,hh4,hh5},
Dh1=Table[Null,{k,1,t}];
Dh2=Table[Null,{k,1,t}];
Dh3=Table[Null,{k,1,t}];
Dh4=Table[Null,{k,1,t}];
Dh5=Table[Null,{k,1,t}];
Dh1[[1]]=0;
Dh2[[1]]=0;
Dh3[[1]]=0;
Dh4[[1]]=0;
Dh5[[1]]=0;
rule={nw2->0,ϵ->0,h1[nw2,ϵ]->0,h2[nw2,ϵ]->0,h3[nw2,ϵ]->0,h4[nw2,ϵ]->0,h5[nw2,ϵ]->0,Derivative[1,0][h1][nw2,ϵ]->0,Derivative[1,0][h2][nw2,ϵ]->0,Derivative[1,0][h3][nw2,ϵ]->0,Derivative[1,0][h4][nw2,ϵ]->0,Derivative[1,0][h5][nw2,ϵ]->0};
For [ict=2,ict≤t, ict++, 
rule=Join[Flatten@Solve[{(Flatten@D[N1[nw2,ϵ]=={0,0,0,0,0},{nw2,ict}]//.rule)},{(D[h1[nw2,ϵ],{nw2,ict}]/.{nw2->0,ϵ->0}),(D[h2[nw2,ϵ],{nw2,ict}]/.{nw2->0,ϵ->0}),(D[h3[nw2,ϵ],{nw2,ict}]/.{nw2->0,ϵ->0}),(D[h4[nw2,ϵ],{nw2,ict}]/.{nw2->0,ϵ->0}),(D[h5[nw2,ϵ],{nw2,ict}]/.{nw2->0,ϵ->0})}],rule];
Dh1[[ict]]=(D[h1[nw2,ϵ],{nw2,ict}]//.rule);
Dh2[[ict]]=(D[h2[nw2,ϵ],{nw2,ict}]//.rule);
Dh3[[ict]]=(D[h2[nw2,ϵ],{nw2,ict}]//.rule);
Dh4[[ict]]=(D[h2[nw2,ϵ],{nw2,ict}]//.rule);
Dh5[[ict]]=(D[h2[nw2,ϵ],{nw2,ict}]//.rule)];

hh1=Sum[1/Factorial[i]*Dh1[[i]]*nw2^{i},{i,1,t}];
hh2=Sum[1/Factorial[i]*Dh2[[i]]*nw2^{i},{i,1,t}];
hh3=Sum[1/Factorial[i]*Dh3[[i]]*nw2^{i},{i,1,t}];
hh4=Sum[1/Factorial[i]*Dh4[[i]]*nw2^{i},{i,1,t}];
hh5=Sum[1/Factorial[i]*Dh5[[i]]*nw2^{i},{i,1,t}];
{nu1->hh1,nu2->hh2,nv1->hh3,nv2->hh4,nw1->hh5}
];
Flatten@Approx[2] //Simplify

Solve[D[N1[nw2,ϵ]=={0,0,0,0,0},{nw2,2}]/.{nw2->0,ϵ->0,h1[nw2,ϵ]->0,h2[nw2,ϵ]->0,h3[nw2,ϵ]->0,h4[nw2,ϵ]->0,h5[nw2,ϵ]->0,Derivative[1,0][h1][nw2,ϵ]->0,Derivative[1,0][h2][nw2,ϵ]->0,Derivative[1,0][h3][nw2,ϵ]->0,Derivative[1,0][h4][nw2,ϵ]->0,Derivative[1,0][h5][nw2,ϵ]->0},{Derivative[2,0][h1][0,0],Derivative[2,0][h2][0,0],Derivative[2,0][h3][0,0],Derivative[2,0][h4][0,0],Derivative[2,0][h5][0,0]}]//Simplify

(* Center manifold approximation *)

NSystem[nu1,nw2,nu2,nv1,nv2,nw1][[2]]/.Approx[2]//Simplify

(* Plot of the pitchfork bifurcation *)

ContourPlot[{Re[{(nw2 (4 z1 ϵ ((ζ1)^2) (Sqrt[-8+ρ1]) ((ρ1)^2) (2+3 ρ1)-(nw2^2) (1+2 ϵ) (-2+ρ1) (-2 (Sqrt[-8+ρ1])-10 (Sqrt[ρ1])+(Sqrt[-8+ρ1]) ρ1+((ρ1)^(3/2)))))/(2 z1 ((ζ1)^2) (Sqrt[-8+ρ1]) ((ρ1)^2) (2+3 ρ1))/.{ζ1->-.5,ζ2->.5,ρ1->.9,ρ2->.9,z1->1,z2->1}}]==0},{ϵ,-.5,.5},{nw2,-1,1},AxesLabel->{ϵ,u}]

(* Pitchfork Bifurcation Conditions - See page 35 of Gibbs (2015) *)
((nw2 (4 z1 ϵ ((ζ1)^2) (Sqrt[-8+ρ1]) ((ρ1)^2) (2+3 ρ1)-(nw2^2) (1+2 ϵ) (-2+ρ1) (-2 (Sqrt[-8+ρ1])-10 (Sqrt[ρ1])+(Sqrt[-8+ρ1]) ρ1+((ρ1)^(3/2)))))/(2 z1 ((ζ1)^2) (Sqrt[-8+ρ1]) ((ρ1)^2) (2+3 ρ1)))/.{nu1->0,nu2->0,nv1->0,nv2->0,nw1->0,nw2->0,ϵ->0}

D[((nw2 (4 z1 ϵ ((ζ1)^2) (Sqrt[-8+ρ1]) ((ρ1)^2) (2+3 ρ1)-(nw2^2) (1+2 ϵ) (-2+ρ1) (-2 (Sqrt[-8+ρ1])-10 (Sqrt[ρ1])+(Sqrt[-8+ρ1]) ρ1+((ρ1)^(3/2)))))/(2 z1 ((ζ1)^2) (Sqrt[-8+ρ1]) ((ρ1)^2) (2+3 ρ1))),nw2]/.{nu1->0,nu2->0,nv1->0,nv2->0,nw1->0,nw2->0,ϵ->0}

D[((nw2 (4 z1 ϵ ((ζ1)^2) (Sqrt[-8+ρ1]) ((ρ1)^2) (2+3 ρ1)-(nw2^2) (1+2 ϵ) (-2+ρ1) (-2 (Sqrt[-8+ρ1])-10 (Sqrt[ρ1])+(Sqrt[-8+ρ1]) ρ1+((ρ1)^(3/2)))))/(2 z1 ((ζ1)^2) (Sqrt[-8+ρ1]) ((ρ1)^2) (2+3 ρ1))),{nw2,2}]/.{nu1->0,nu2->0,nv1->0,nv2->0,nw1->0,nw2->0,ϵ->0}

D[((nw2 (4 z1 ϵ ((ζ1)^2) (Sqrt[-8+ρ1]) ((ρ1)^2) (2+3 ρ1)-(nw2^2) (1+2 ϵ) (-2+ρ1) (-2 (Sqrt[-8+ρ1])-10 (Sqrt[ρ1])+(Sqrt[-8+ρ1]) ρ1+((ρ1)^(3/2)))))/(2 z1 ((ζ1)^2) (Sqrt[-8+ρ1]) ((ρ1)^2) (2+3 ρ1))),{nw2,3}]/.{nu1->0,nu2->0,nv1->0,nv2->0,nw1->0,nw2->0,ϵ->0}

D[((nw2 (4 z1 ϵ ((ζ1)^2) (Sqrt[-8+ρ1]) ((ρ1)^2) (2+3 ρ1)-(nw2^2) (1+2 ϵ) (-2+ρ1) (-2 (Sqrt[-8+ρ1])-10 (Sqrt[ρ1])+(Sqrt[-8+ρ1]) ρ1+((ρ1)^(3/2)))))/(2 z1 ((ζ1)^2) (Sqrt[-8+ρ1]) ((ρ1)^2) (2+3 ρ1))),{ϵ}]/.{nu1->0,nu2->0,nv1->0,nv2->0,nw1->0,nw2->0,ϵ->0}

D[((nw2 (4 z1 ϵ ((ζ1)^2) (Sqrt[-8+ρ1]) ((ρ1)^2) (2+3 ρ1)-(nw2^2) (1+2 ϵ) (-2+ρ1) (-2 (Sqrt[-8+ρ1])-10 (Sqrt[ρ1])+(Sqrt[-8+ρ1]) ρ1+((ρ1)^(3/2)))))/(2 z1 ((ζ1)^2) (Sqrt[-8+ρ1]) ((ρ1)^2) (2+3 ρ1))),{{nw2,ϵ},2}]/.{nu1->0,nu2->0,nv1->0,nv2->0,nw1->0,nw2->0,ϵ->0}