(* Proof of Lemma 1 Gibbs (2015); Last Updated 25/09/2015 *)
(* The Tmap *)
T[M1_List,M2_List,w_List]:=Module[{A,B,CC,R,G,Mhold,Ma,Mb,TM1,TM2,TW},
A=μ+α (w[[1]]M1[[1]]+w[[2]]M2[[1]]);
B=α w[[1]]M1[[2]]+ζ1;
CC=α w[[2]] M2[[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}};
R={{A},{M1[[1]]},{M2[[1]]}};
G={{B,CC},{M1[[2]],0},{0,M2[[2]]}};
Mhold=R.R\[Transpose]+G.{{z1,s},{s,z2}}.G\[Transpose];

TW={{(Mhold[[1]][[2]]-w[[2]] Mhold[[2]][[3]])/Mhold[[2]][[2]]},{(Mhold[[1]][[3]]-w[[1]]Mhold[[2]][[3]])/Mhold[[3]][[3]]}};
{TM1[[1]]-M1[[1]],TM1[[2]]-M1[[2]],TM2[[1]]-M2[[1]],TM2[[2]]-M2[[2]],TW[[1]]-w[[1]],TW[[2]]-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=b1 (-1+α γ1)+ζ1;
ss3=b2 (-1+α γ2)+ζ2;
ss4=(1/((a1^2)+(b1^2) z1))((a1^2) (-1+α) γ1+a1 a2 (-1+α) γ2+b1 z1 (b1 (-1+α) γ1+ζ1));
ss5=(1/((a2^2)+(b2^2) z2))(a1 a2 (-1+α) γ1+(a2^2) (-1+α) γ2+b2 z2 (b2 (-1+α) γ2+ζ2));
{cc,ss1,ss2,ss3,ss4,ss5}
];
TransFormSys[u1,u2,v1+ζ1/(1-(ϵ+1/2)),v2+ζ2/(1-(ϵ+1/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)),v2+ζ2/(1-(ϵ+1/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]]

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

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)/2)+v1 ϵ+v1 w1 ϵ-((ζ1)/(2 ((1/2)-ϵ)))+((w1 ζ1)/(2 ((1/2)-ϵ)))+((ϵ ζ1)/((1/2)-ϵ))+((w1 ϵ ζ1)/((1/2)-ϵ));
f4=((v2 w2)/2)+v2 ϵ+v2 w2 ϵ-((ζ2)/(2 ((1/2)-ϵ)))+((w2 ζ2)/(2 ((1/2)-ϵ)))+((ϵ ζ2)/((1/2)-ϵ))+((w2 ϵ ζ2)/((1/2)-ϵ));
f5=((u1^2) (1+w1) (-(1/2)+ϵ)+u1 u2 (1+w2) (-(1/2)+ϵ)+1/(-2+4 ϵ)z1 (v1 (-1+2 ϵ)-2 ζ1) (v1 (1+w1) (-1+2 ϵ)-2 w1 ζ1))/((u1^2)+z1 (((v1+((2 ζ1)/(1-2 ϵ))))^2));
f6=(u1 u2 (1+w1) (-(1/2)+ϵ)+(u2^2) (1+w2) (-(1/2)+ϵ)+1/(-2+4 ϵ)z2 (v2 (-1+2 ϵ)-2 ζ2) (v2 (1+w2) (-1+2 ϵ)-2 w2 ζ2))/((u2^2)+z2 (((v2+((2 ζ2)/(1-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+w2,u1+w2,2 I v1 ζ1-2 I w1 ζ1,2 I u2 ζ2-2 I v2 ζ2,v1+w1,u2+v2]

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

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

(* 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]

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,nu2,nv1,nv2,nw1,nw2][[6]]/.Approx[2]

(* Plot of the pitchfork bifurcation *)

ContourPlot[{Re[{2 nw2 ϵ-((((1/12)-(I/6)) (nw2^3))/(z2 ((ζ2)^2)))-((((1/6)-(I/3)) (nw2^3) ϵ)/(z2 ((ζ2)^2))) /.{ζ1->.5,ζ2->.5,z1->1,z2->1}}]==0},{ϵ,-.5,.5},{nw2,-1,1},AxesLabel->{ϵ,u}]

(* Pitch Fork Bifurcation Conditions. See Page 35 of Gibbs (2015) *)
2 nw2 ϵ-((((1/12)-(I/6)) (nw2^3))/(z2 ((ζ2)^2)))-((((1/6)-(I/3)) (nw2^3) ϵ)/(z2 ((ζ2)^2)))/.{nu1->0,nu2->0,nv1->0,nv2->0,nw1->0,nw2->0,ϵ->0}

D[2 nw2 ϵ-((((1/12)-(I/6)) (nw2^3))/(z2 ((ζ2)^2)))-((((1/6)-(I/3)) (nw2^3) ϵ)/(z2 ((ζ2)^2))),nw2]/.{nu1->0,nu2->0,nv1->0,nv2->0,nw1->0,nw2->0,ϵ->0}

D[2 nw2 ϵ-((((1/12)-(I/6)) (nw2^3))/(z2 ((ζ2)^2)))-((((1/6)-(I/3)) (nw2^3) ϵ)/(z2 ((ζ2)^2))),{nw2,2}]/.{nu1->0,nu2->0,nv1->0,nv2->0,nw1->0,nw2->0,ϵ->0}

D[2 nw2 ϵ-((((1/12)-(I/6)) (nw2^3))/(z2 ((ζ2)^2)))-((((1/6)-(I/3)) (nw2^3) ϵ)/(z2 ((ζ2)^2))),{nw2,3}]/.{nu1->0,nu2->0,nv1->0,nv2->0,nw1->0,nw2->0,ϵ->0}

D[2 nw2 ϵ-((((1/12)-(I/6)) (nw2^3))/(z2 ((ζ2)^2)))-((((1/6)-(I/3)) (nw2^3) ϵ)/(z2 ((ζ2)^2))),{ϵ}]/.{nu1->0,nu2->0,nv1->0,nv2->0,nw1->0,nw2->0,ϵ->0}

D[2 nw2 ϵ-((((1/12)-(I/6)) (nw2^3))/(z2 ((ζ2)^2)))-((((1/6)-(I/3)) (nw2^3) ϵ)/(z2 ((ζ2)^2))),{{nw2,ϵ},2}]/.{nu1->0,nu2->0,nv1->0,nv2->0,nw1->0,nw2->0,ϵ->0}