# To solve an ode system with pulses

I want to solve an ode system with pulses, for example

x'[t] == b*(1 - x[t]) - a*x[t]*y[t],
y'[t] == a*x[t]*y[t] - d*y[t]
in {t,0,1}


and then at the end of each period, a pulse occurs, for example, in the next period in {t,1,2},

x[0]=x[1],y[0]=y[1]*(1+x[1]),...,


and so on. I made the following code but it seems not to work.

IDESolver[{x0_, y0_}, {a_, d_, b_, T_},
MaxPeriod_] := Module[{sol, resulttable},
resulttable = {};
Do[
sol =
NDSolve[{x'[t] == b*(1 - x[t]) - a*x[t]*y[t],
y'[t] == a*x[t]*y[t] - d*y[t], x[0] == x0,
y[0] == y0}, {x, y}, {t, 0, T}];
AppendTo[resulttable, sol];
x[0] = Evaluate[x[T] /. sol];
y[0] = Evaluate[(y[T]*(1+x[T])) /. sol] ;, {i, 1, MaxPeriod}
]
resulttable
]

IDESolver[{0.2, 0.2}, {0.3, 0.1, 0.2, 1}, 10]


Can anybody help me with this issue? Thanks in advance.

• "It seems not to work..." What about it doesn't work? Does it spit out errors? Does it spit out a wrong solution (if so, what is the solution you expect and how does it differ from the wrong solution)? Have you narrowed down where the problem is? Etc. – march Dec 18 '15 at 16:59
• It does not work and gives errors such as "$RecursionLimit::reclim: Recursion depth of 256 exceeded. >>" – X Bruno Dec 18 '15 at 17:03 ## 2 Answers I made a small change and i think it is now working.  IDESolver[{x0_, y0_}, {a_, d_, b_, T_}, MaxPeriod_] := Block[{X0 = x0, Y0 = y0, sol, resulttable}, resulttable = {}; Do[ sol = NDSolve[{x'[t] == b*(1 - x[t]) - a*x[t]*y[t], y'[t] == a*x[t]*y[t] - d*y[t], x[0] == X0, y[0] == Y0}, {x, y}, {t, 0, T}]; AppendTo[resulttable, sol]; X0 = Evaluate[x[T] /. sol]; Y0 = Evaluate[(y[T]*(1 + x[T])) /. sol];, {i, 1, MaxPeriod}] ; Flatten[resulttable] ]; listsols = IDESolver[{0.2, 0.2}, {0.3, 0.1, 0.2, 1}, 10]; Table[Plot[listsols[[i]][[2]][x], {x, 0, 10}], {i, 1, Length[listsols]}]  • Thank you very much! It still gives errors such as "$RecursionLimit::reclim: ""Recursion depth of 256 exceeded." PS. I am using Mathematica 7 on Mac. – X Bruno Dec 18 '15 at 18:05
• I am using Mathematica 10 on Windows 8 and its is working fine... – Diogo Dec 18 '15 at 18:08
• I restart the Kernel and it is working! Thanks a lot! – X Bruno Dec 18 '15 at 18:16
• You are welcome! – Diogo Dec 18 '15 at 18:18

In the case where you actually want a continuous time variable, the function WhenEvent will be extremely useful. I've modified your code to:

Clear[IDESolver, x, y, sol]
IDESolver[{x0_, y0_}, {a_, d_, b_, period_}, tMax_] :=
NDSolve[{
x'[t] == b*(1 - x[t]) - a*x[t]*y[t]
, y'[t] == a*x[t]*y[t] - d*y[t]
, x[0] == x0, y[0] == y0
, WhenEvent[Mod[t, period] == 0, y[t] -> y[t] (1 + x[t])]
}
, {x, y}
, {t, 0, tMax}]


Then,

sol = First@IDESolver[{0.2, 0.2}, {0.3, 0.1, 0.2}, 1, 10];
Plot[Evaluate[# /. sol], {t, 0, 10}] & /@ {x[t], y[t]} // GraphicsRow


• @SteveX. You're welcome! Don't forget to accept the answer that you find most helpful by clicking the grey checkmark next to that answer! – march Dec 21 '15 at 19:23