Using odeint to determine bifurcation

Question

I am attempting to get the data for the bifurcation diagram, namely the values of trigger [A] (which are specified by linspace of Aspace) which give steady-state (nullcline) values of [N] in this set of differential equations:

As the set of equations is too long to attempt to get the nullclines by re-arranging, I am trying to use odeint to integrate over given time, and get the values of all d[x]/d(t) to 0, and return the values of [A] and [N] for those particular nullclines. In theory, this code should return the values expected, however it seems to kill the kernel without further explanation, preventing me from tracing down the issue. Please find the code below (don't get scared, most of it is just stating the variable values).

My code:

#basic conditions as specified in the material
OS = 0
O = 0
S= 0
N = 0
B = 0

n1 = 1*(10**-4)
a1 = 1
a2 = 0.01
a3 = 0.2
n2 = 1*(10**-7)
b1 = 0.0011
b2 = 0.001
b3 = 0.0007
n3 = 1*(10**-4)
c1 = 1
c2 = 0.01
c3 = 0.2
n4 = 1*(10**-7)
d1 = 0.0011
d2 = 0.001
d3 = 0.0007
n5 = 1*(10**-4)
e1 = 0.005
e2 = 0.1
n6 = 1*(10**-7)
f1 =0.001
f2 = 9.95*(10**-4)
f3 = 0.01
k1c = 0.05
k2c = 0.001
k3c = 5
g1 = 1
g2 = 1
g3 = 1
def differential_eq(y, t, A, B):
    O, S, OS, N = y
    dydt = np.empty(4)
    dydt[0] = ((n1 + a1*A + a2*OS + a3*OS*N)/(1 + n2 + b1*A + b2*OS + b3*OS*N)) - g1*O - k1c*O*S + k2c*OS
    dydt[1] = (n3 + c1*A + c2*OS + c3*OS*N)/(1 + n4 + d1*A + d2*OS + d3*OS*N) - g2*S - k1c*O*S + k2c*OS
    dydt[2] = k1c*O*S -k2c*OS - k3c*OS
    dydt[3] = ((n5 + e1*OS + e2*OS*N)/(1 + n6 +f1*OS + f2*OS*N + f3*B)) - g3*N
    return dydt

timepool =np.linspace(0,100,10)

def simulate(A):
A = A
y = (0,0,0,0)
B =0
solution = sp.odeint (differential_eq(y, timepool, A, B), (A, B), timepool)
    solution = sp.odeint (differential_eq, initial, timepool)
    if dydt[0] == 0 and dydt[1] ==0 and dydt[2] ==0 and dydt[3] ==0:
        print (A,N, solution)
        return (A, N)

Aspace = np.linspace(0,150,150)
for A in Aspace:
    simulate(A)

Now, this does not seem to work, and I see no indication of why, as it just kills the kernel: if anyone has any idea of why this is the case, please let me know.

Answer 1

Combining the remarks in the comments and some more, your simulate function should look like

def simulate(A):
    y = (0,0,0,0)
    B =0
    # parameter for the accuracy level
    eps = 1e-9
    # pass the correct arguments: first the function reference, then
    # initial point and sample times, the additional arguments and
    # the absolute and relative tolerances (they are combined as atol+rtol*(norm(y))
    solution = sp.odeint (differential_eq,y, timepool, args = (A, B), atol=1e-4*eps, rtol=1e-3*eps)
    y, t = solution[-1], timepool[-1]
    # compute the slopes at the last sample point by calling the ODE function
    dydt = differential_eq(y,t,A,B)
    # variables not explicitly declared global are local to the ODE function, 
    # they do not influence the global variables. 
    O, S, OS, N = y
    # never compare floating point numbers with "==" in numerics
    # always account for the floating point noise accumulated in the computation
    if max(abs(dydt)) < eps :
        #here print only the last sample point 
        print "%6.2f, %16.12f, %s"%(A, N, y)

    return (A, N)

which gives the result (shortened)

  0.00,   0.000099999991, [  9.99994911e-05   9.99994911e-05   9.99789864e-11   9.99999905e-05]
 10.00,   0.002883872464, [  7.25813895e+00   7.25813895e+00   5.26700470e-01   2.88387246e-03]
 20.00,   0.008780897525, [  1.21628435e+01   1.21628435e+01   1.47905181e+00   8.78089753e-03]
 30.00,   0.017500871488, [ 16.0785845   16.0785845    2.58469186   0.01750087]
 40.00,   0.030166487759, [ 19.40580443  19.40580443   3.76509943   0.03016649]
 50.00,   0.049390699930, [ 22.32996985  22.32996985   4.98527848   0.0493907 ]
 60.00,   0.081339098240, [ 24.95672023  24.95672023   6.22713342   0.0813391 ]
 70.00,   0.144113671629, [ 27.35699725  27.35699725   7.48255648   0.14411367]
100.00,  63.103453906658, [ 52.49789017  52.49789017  27.55477377  63.10345391]
110.00,  64.363708289042, [ 53.43024046  53.43024046  28.54219752  64.36370829]
120.00,  65.425943398185, [ 54.25586731  54.25586731  29.43110515  65.4259434 ]
130.00,  66.343605498371, [ 55.00078061  55.00078061  30.24480971  66.3436055 ]
140.00,  67.150809538919, [ 55.68201657  55.68201657  30.99866996  67.15080954]
150.00,  67.870747329663, [ 56.31143752  56.31143752  31.70343926  67.87074733]