Author 
Comment/Response 
Derek

06/28/08 12:05pm
Hello,
I've been trying to generate bifurcation diagrams for a system of differential equations but got stuck in the technicalities.
If anyone can show me how to do it simply then i'd be very grateful. I post my attempt below (from copying an example from a simpler system). if someone could show me where i'm going wrong then that would be great!! I would like to vary a parameter, for example v, and see how the dynamics change over time is this doable? I have seen it done for discrete single species systems but cannot find anything for something like this. Many thanks, Derek
system::::==================================
[{If[H[t] < 0, 0, a H[t]  ap P[t] + at T[t]  (H[t] b + bp P[t] + bt T[t])] == H'[t], If[P[t] < 0  H[t] < 0, 0, P[t] ((L H[t])/(H0 + H[t])  u)  (H[t] b + bp P[t] + bt T[t]) (P[t]/H[t]) + phi T[t]] == P'[t], If[T[t] < 0  H[t] < 0, 0,
H[t] (theta  rho)  (H[t] b + bp P[t] + bt T[t]) (T[t]/H[t])  v P[t]] == T'[t],
==============================================
attempt:::
a = 3; b = 1.31; u = .1; H0 = 10; ap = .3; L = 4; bp = 1.2;
at = .2; bt = 1.2; phi = .8; theta = .75; rho = .005; v = .213; tstop
= 50;
tmin = 0; tmax = 50; H0 = 5; P0 = 5; T0 = 5;
sol[v] = NDSolve[{If[H[t] < 0, 0,
a H[t]  ap P[t] + at T[t]  (H[t] b + bp P[t] + bt T[t])] ==
H'[t], If[P[t] < 0  H[t] < 0, 0,
P[t] ((L H[t])/(H0 + H[t])  u)  (H[t] b + bp P[t] +
bt T[t]) (P[t]/H[t]) + phi T[t]] == P'[t],
If[T[t] < 0  H[t] < 0, 0,
H[t] (theta  rho)  (H[t] b + bp P[t] + bt T[t]) (T[t]/H[t]) 
v P[t]] == T'[t], H[tmin] == H0, P[tmin] == P0,
T[tmin] == T0}, {H[t], T[t], P[t]}, {t, tmin, tmax}];
Clear[bifurcationplot]
bifurcationplot[{vbegin_, vend_, vstep_}, timesteps_, takelast_] :=
Module[{H, iter, v, thelist = {}}, H[v_] := Evaluate[H[t] /. sol[v]]
iter[v_] := Table[sol[v], {t, 0, timesteps}];
For[v = vbegin, v < vend, v += vstep,
thelist =
Flatten[{thelist, ({v, #1} &) /@ Take[iter[v], takelast]}, 1]];
ListPlot[thelist, PlotStyle > PointSize[0.002],
PlotRange > {Automatic, {0, 2000}}]]
bifurcationplot[{0, 10, 0.3}, 200, 20]
URL: , 
