Why is the integration not progressing with scipy.integrate.ode?

Question

I am integrating a function func in the complex domain using scipy.integrate.ode. I have assembled the following code structure:

import numpy as np
from scipy.integrate import ode

def func(x, u, k):
    #f, g = u
    dfdx = k*x**2
    dgdx = -x
    
    rhs_FD = [0]*2
    rhs_FD[0] = dfdx
    rhs_FD[1] = dgdx

    return np.array(rhs_FD)

solver = ode(func)
solver.set_integrator('zvode',method='bdf')

k = 1+1j
solver.set_f_params(k)

x0 = 0.0
u0 = [0.0, 1.0]
solver.set_initial_value(u0, x0)

x_values = np.linspace(x0, 1.0, 3)

u_values = []

for x in x_values:
    u = solver.integrate(x)
    u_values.append(u)

However, I noticed that the choice of methods, described by solver.set_integrator('zvode', method='bdf'), is not allowing the integration to advance. The output remains the same as the initial condition:

[[0.+0.j 1.+0.j]
 [0.+0.j 1.+0.j]
 [0.+0.j 1.+0.j]]

I noticed that the code works well by replacing the integrator with solver.set_integrator('dopri5') and considering a real value of k. However, I need to perform this integration in the complex domain. How can I resolve the issue of integration not advancing when considering the complex domain?

Answer 1

If all you care about is getting this working then the issue is that you already initialized the solver but you are asking it to evaluate it again at x=0. To fix this, loop through x_values[1:].

for x in x_values[1:]:
    u = solver.integrate(x)
    u_values.append(u)
    
print(u_values)

Output:

[array([0.04166667+0.04166667j, 0.87499999+0.j        ]), array([0.33333333+0.33333333j, 0.5       +0.j        ])]