PFLocation[alpha_,eta_] := (

If[ 1 + alpha eta < 0.0,
    (* the pitchfork is on the in-phase branch *)
    gammaPF = -(1.0+alpha eta)/(1.0+alpha^2),
    (* the pitchfork is on the anti-phase branch *)
    gammaPF = (1.0+alpha eta)/(3.0+alpha^2+2.0 alpha eta) ];
    
{eta, gammaPF} )

PlotPFLocation[alpha_,etamin_,etamax_,etastep_] := (

temp = Table[ PFLocation[alpha,eta], 
              {eta,etamin,etamax,etastep} ];
ListPlot[temp] )

PlotPFLocation[-2,-5,5,0.1]

PFCriticality[alpha_,eta_] := (

(* For the given values of alpha and eta,
   determine the criticality of the pitchfork bifurcation 
   occurring on either the inphase or antiphase branch of
   steady states of the reduced system (for identical cells). *)

(* The reduced system at Delta=1 *)
h1[r1_,r2_,phi_,gamma_] := r1 (1-r1^2) + 
                           gamma (Cos[phi] - alpha Sin[phi]) r2 -
                           gamma r1;
h2[r1_,r2_,phi_,gamma_] := r2 (1-r2^2) +
                           gamma (Cos[phi] + alpha Sin[phi]) r1 - 
                           gamma r2;
h3[r1_,r2_,phi_,gamma_] := eta (r1^2-r2^2) -
                           gamma r2/r1 (Sin[phi] + alpha Cos[phi]) - 
                           gamma r1/r2 (Sin[phi] - alpha Cos[phi]);
                           
If[ 1 + alpha eta < 0.0,
    (* the pitchfork is on the in-phase branch *)
    gammaPF = -(1.0+alpha eta)/(1.0+alpha^2); 
    r10 = 1.0; 
    r20 = 1.0; 
    phi0 = 0.0,
    (* the pitchfork is on the anti-phase branch *)
    gammaPF = (1.0+alpha eta)/(3.0+alpha^2+2.0 alpha eta);
    r10 = Sqrt[1.0-2.0 gammaPF];
    r20 = Sqrt[1.0-2.0 gammaPF];
    phi0 = N[Pi,10] ]; 
    
vars = {r1, r2, phi};
matD[F_,X_] := Outer[D,F,X];

jac = matD[ {h1[r1,r2,phi,gamma],
             h2[r1,r2,phi,gamma],
             h3[r1,r2,phi,gamma]}, vars ];
l0 = jac /. {r1->r10, r2->r20, phi->phi0, gamma->gammaPF};
l1 = D[jac,gamma] /. {r1->r10, r2->r20, phi->phi0, gamma->gammaPF};
m = matD[jac,vars]/2;
m0 = m /. {r1->r10, r2->r20, phi->phi0, gamma->gammaPF};
n = matD[m,vars]/3;
n0 = n /. {r1->r10, r2->r20, phi->phi0, gamma->gammaPF};

{vals,vecs} = Eigensystem[l0];
For[ i=1, i<=3 && Abs[vals[[i]]] > 0.0001, i++ ];
v0 = vecs[[i]];
indices = Complement[{1,2,3},{i}];
v1 = vecs[[Part[indices,1]]];
v2 = vecs[[Part[indices,2]]];
v0 /= Sqrt[v0.v0];
v1 /= Sqrt[v1.v1];
v2 /= Sqrt[v2.v2];
evalue1 = vals[[Part[indices,1]]];
evalue2 = vals[[Part[indices,2]]];

{vals,vecs} = Eigensystem[Transpose[l0]];
For[ i=1, i<=3 && Abs[vals[[i]]] > 0.0001, i++ ];
v0star = vecs[[i]];
indices = Complement[{1,2,3},{i}];
If[ vals[[Part[indices,1]]] == evalue1,
    v1star = vecs[[Part[indices,1]]];
    v2star = vecs[[Part[indices,2]]],
    v1star = vecs[[Part[indices,2]]];
    v2star = vecs[[Part[indices,1]]] ];
v0star /= v0star.v0;
v1star /= v1star.v1;
v2star /= v2star.v2;

u1 = v0;
gamma1 = -m0.u1.u1.v0star / l1.u1.v0star;

b = -m0.u1.u1;
u2 = v1star.b/evalue1 v1 + v2star.b/evalue2 v2;
gamma2 = -(2 m0.u1.u2 + n0.u1.u1.u1).v0star / (v0star.l1.u1);

(* Return gamma2.  
   If gamma2 > 0 (gamma2 < 0), the pitchfork bifurcation is 
   subcritical (supercritical). *)
{eta,gamma2} )

PlotPFCriticality[alpha_,etamin_,etamax_,etastep_] := (

etanum = Ceiling[ (etamax-etamin)/etastep ];
temp = Table[ PFCriticality[alpha,etamin+etastep i], 
              {i,0,etanum} ];
ListPlot[temp] )

NeutralCrit[alpha_,etamin_,etamax_,tol_] := (

(* For the given value of alpha, use bisection to 
   find the value of eta (within tolerance tol) at 
   which the pitchfork bifurcation has neutral criticality. *)
   
fmid = Part[PFCriticality[alpha,etamax], 2];
f    = Part[PFCriticality[alpha,etamin], 2];
If[ f fmid >= 0, 
    (* Root must be bracketed *)
    {alpha,Null},
    (* Go ahead with the bisection *)
    If[ f <= 0, 
        etabis = etamin; deta = etamax - etamin,
        etabis = etamax; deta = etamin - etamax ];
    While[ !(Abs[deta] < tol || fmid==0),
           deta = deta / 2;
           etamid = etabis + deta;
           fmid = Part[PFCriticality[alpha,etamid],2];
           If[ fmid <= 0, etabis = etamid ] ];
    N[{alpha,etabis}] ] )

PlotNeutralCrit[alphamin_,alphamax_,alphastep_,
                mineta_,maxeta_,tol_] := (
             
(* *)

temp = {};
Do[ (* If[ region==1,
        etamin = -1/alpha+0.05; etamax = maxeta,
        If[ region==2,
        etamin = (alpha^2-1)/(2 alpha); etamax = (-1/alpha-0.05),
        etamin = mineta; etamax = (alpha^2-1)/(2 alpha) ] ]; *)
    etamin = mineta;
    etamax = maxeta;
    temp = Append[temp, NeutralCrit[alpha,etamin,etamax,tol] ],
    
    {alpha,alphamin,alphamax,alphastep} ];
ListPlot[temp] )

PrintTable[table_] := (
out = OpenWrite["neutralpf.out"];
Do[ WriteString[out, ToString[Part[Part[table,i],1]],
                     "  ",
                     ToString[Part[Part[table,i],2]] 
                     "\n"], 
    {i,1,Length[table]} ];
Close[out]; )

PFCriticality[-0.522,15.555]

PlotPFCriticality[-1.4,-5,5,0.2]

NeutralCrit[1.65,4,6,0.01]

PlotNeutralCrit[0.05,2.05,0.1,-6,-1.2,0.01]

PrintTable[temp]