(*
 * predprey_bif.mma - predator prey bifurcation analysis
 *
 * RMM, 8 Jul 06
 *)

<<Jac.m
<<predprey.mma

nompar = {r->1.6, d->0.56, b->0.6, k->125, a->3.2, c->50}
bifpar = {r->1.6, d->0.56, b->0.6, k->125}

App = FullSimplify[ Jac[Fpp[H, L], {H, L}] /. eqpts[[1]] ];
eigpp = FullSimplify[ Eigenvalues[App] ];

(* Create a function to compute the real part of the largest eigenvalue *)
maxeig[aval_, cval_] :=
  Max @@ Thread[Re[eigpp /. {a->aval, c->cval} /. bifpar]];

(* Compute equilibrium points/stability as a function of a *)
eqval1 = {}; eqval2 = {}; eqval1 = {};
stab1 = {}; stab2 = {}; stab1 = {};
avec = {}; mxval1 = {}; mnval1 = {};
cmin={}; cmax={}; 
For[aval = 1.35; ind=1, aval < 10, aval += 0.05; ind++,
  (* Compute the equilibrium points for these parameter values *)
  AppendTo[eqval1, H /. eqpts[[1]] /. a->aval /. nompar];
  AppendTo[eqval2, H /. eqpts[[2]] /. a->aval /. nompar];
  AppendTo[eqval3, H /. eqpts[[3]] /. a->aval /. nompar];

  (* Determine the stability of the equilibrium point *)
  AppendTo[stab1, Max @@ Thread[Re[eigpp /. a->aval /. nompar]]];

  (* Simulate the system and determine excursions *)
  L0 = L /. eqpts[[1]] /. a->aval /. nompar;
  H0 = H /. eqpts[[1]] /. a->aval /. nompar;
  soln = NDSolve[
    { H'[t] == Fpp[H[t], L[t]][[1]], L'[t] == Fpp[H[t], L[t]][[2]],
      H[0] == H0 + 5, L[0] == L0 + 5 } /. a->aval /. nompar,
    {H, L}, {t, 0, 250}
  ];
  tvec = Array[#&, 500]/10 + 200;
  AppendTo[mxval1, Max @@ (MapThread[H, {tvec}] /. soln)];
  AppendTo[mnval1, Min @@ (MapThread[H, {tvec}] /. soln)];

  (* Figure out some stability boundaries *)
  AppendTo[cmin, c /. FindRoot[maxeig[aval, c] == 0, {c, 10, 10, 110}]];
  AppendTo[cmax, c /. FindRoot[maxeig[aval, c] == 0, {c, 150, 50, 200}]];

  AppendTo[avec, aval];
];

(* Generate a bifurcation diagram *)
eqplot = ListPlot[Transpose[{avec, eqval1}], PlotJoined->True, PlotStyle->Thickness[0.01]];
mxplot = ListPlot[Transpose[{avec, mxval1}], PlotJoined->True, PlotStyle->Dashing[{0.025, 0.025}]];
mnplot = ListPlot[Transpose[{avec, mnval1}], PlotJoined->True, PlotStyle->Dashing[{0.025, 0.025}]];
Show[{eqplot, mxplot, mnplot}];

(* Generate simulations corresponding to a couple of different parameters *)
L0 = L /. eqpts[[1]] /. a->1 /. nompar;
H0 = H /. eqpts[[1]] /. a->1 /. nompar;
soln2 = NDSolve[
  { H'[t] == Fpp[H[t], L[t]][[1]], L'[t] == Fpp[H[t], L[t]][[2]],
    H[0] == H0 + 0.5, L[0] == L0 - 0.5 } /. a->1 /. nompar,
  {H, L}, {t, 0, 100}
];
par2 = Plot[{H[t] /. soln2, L[t] /. soln2}, {t, 0, 100}];

L0 = L /. eqpts[[1]] /. a->3 /. nompar;
H0 = H /. eqpts[[1]] /. a->3 /. nompar;
soln3 = NDSolve[
  { H'[t] == Fpp[H[t], L[t]][[1]], L'[t] == Fpp[H[t], L[t]][[2]],
    H[0] == H0 + 0.5, L[0] == L0 + -0.5 } /. a->3 /. nompar,
  {H, L}, {t, 0, 100}
];
par3 = Plot[{H[t] /. soln3, L[t] /. soln3}, {t, 0, 100}];

L0 = L /. eqpts[[1]] /. a->7.5 /. nompar;
H0 = H /. eqpts[[1]] /. a->7.5 /. nompar;
soln1 = NDSolve[
  { H'[t] == Fpp[H[t], L[t]][[1]], L'[t] == Fpp[H[t], L[t]][[2]],
    H[0] == H0 + 0.5, L[0] == L0 - 0.5 } /. a->7.5 /. nompar,
  {H, L}, {t, 0, 500}
];
par1 = Plot[{H[t] /. soln1, L[t] /. soln1}, {t, 0, 500}];

(* Plot stable versus unstable region *)
ContourPlot[maxeig[x, y], {x, 1.35, 10}, {y, 0, 200}, Contours->{0}, PlotPoints->100];

(* Create a file with the relevant data as well *)
Export["predprey.dat",
  Transpose[{avec, eqval1, stab1, mnval1, mxval1, cmin, cmax}]];