Time-dependent events in ODE

Specific Domains Modelling & Simulations
15

It was fine, you just needed to reduce the tolerance.

using Plots
using DifferentialEquations

i = 1.25 # 0.089;4/45/4=1.25
f = .25
u = .4
w = 1.14
a = 1.19
d = .25
s = .02
b = 0.55
e = .15
r = 5
m = 0
q = .00001
h =1
v = .03

maxtime = 30
tspan = (0.0, maxtime)

stim = 5000

x0 = [10000 1 20.0 200.0 200.0]

function system(dy, y, p, t)

  V = y[1]
  X = y[2]
  Y = y[3]
  P = y[4]
  Q = y[5]

  V, X, Y, P, Q, = y
  T = X+Y
  E = (V)./(v+T+V)
  A = (V)./(s+V)
  B = (V)./(r+T+V)
  C = X./(X+Y+q)
  D = Y./(X+Y+q)

  dy[1] = i - f*V - u*X*A- w*Y*B
  dy[2] = a*X*A - a*X*(1-A) + m*P - m*P*C + h*P*E
  dy[3] = b*Y*B -e*Y - b*Y*(1-B) - m*P*D
  dy[4] = a*X*(1-A) - m*P + m*P*C - d*P - h*P*E
  dy[5] = b*Y*(1-B) + m*P*D - e*Q
end

function condition1(y, t, integrator)
  t - 6
end
function condition2(y, t, integrator)
  t-20
end

function affect!(integrator)
  integrator.u[1] = stim
end

cb1 = ContinuousCallback(condition1, affect!)
cb2 = ContinuousCallback(condition2, affect!)
cb = CallbackSet(cb1,cb2)

prob = ODEProblem(system, x0, tspan)
sol = solve(prob, Rodas4(), callback = cb, reltol=1e-6, isoutofdomain=(u,p,t)->any(x->x<0,u))

plot(sol,layout=(2,3))

savefig("result.png")

result

It was just hard to see what it was doing because of the divergence when values hit negative. You may want to do something like isoutofdomain=(u,p,t)->any(x->x<0,u), so

sol = solve(prob, Rodas4(), callback = cb, isoutofdomain=(u,p,t)->any(x->x<0,u))

also works just fine.

Are you looking for plot(sol,vars=1)? The docs are here: http://docs.juliadiffeq.org/latest/basics/plot.html