integrate.solve_ivp bugged

Question

I am trying to solve an ODE with solve_ivp, but I am getting strange errors. Documentation: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

from scipy.integrate import solve_ivp as solve
import numpy as np
import matplotlib.pyplot as plt

def fun(t, y, beta, eps):
    u = 1/eps * (-1/3 * y[0]**3 + y[0] + y[1])
    v = beta - y[0]
    y[0]= u
    y[1]= v

    return y

eps = 0.05
beta = 1.1
y0 = np.array([2, 2])
t = np.linspace(0, 1000, num=10000)
Y = solve(fun=lambda t, y: fun(t, y, beta, eps), t_span=(t[0],t[-1]), y0=y0 , method='RK45', t_eval=[1,2,3,4]).y

ValueError: need at least one array to concatenate

Note that I already had to deviate from the documentation, as args argument does not work either. Instead, I had to redefine my function which gets 2 additional arguments using the lambda method. I found this here: https://stackoverflow.com/questions/48245765/pass-args-for-solve-ivp-new-scipy-ode-api

Am I doing something wrong or is this method extremely buggy?

Answer 1

Your issue is that you are modifying y in your function. If you return, [u,v] instead of resetting the values in y, it seems to run just fine. I believe this is because the solver will store and use intermediate values elsewhere, but your function is having the side effect of altering these stored values.