Python solve_ivp [scipy.integrate]: how to set the number of output points of your integration process?

Question

I want to integrate the same differential equation with different initial conditions and I want to compare each points that I have obtained. Unfortunately if I do:

sol = solve_ivp(my_model, [0, leng], [1,1,1], method="LSODA", args=(rho, sigma, beta), dense_output=False)

sol_2 = solve_ivp(my_model, [0, leng], [1+0.001,1+0.001, 1-0.001], method="LSODA", args=(rho, sigma, beta), dense_output=False)

The number of points obtained in these two integration processes is different.

I don't want to use the interpolation procedure since I would like to work with real data.

Can I set the number of points of the solution(s)?

Answer 1

One variant to compare instances differing in initial conditions is to collate them all into one big system where they are all solved simultaneously all with the same step sequence

def multi_Lorenz(t,u):
    return np.concatenate([Lorenz(t,uu) for uu in u.reshape([-1,3])])

u0 = np.concatenate([[1+k*1e-5,1+k*1e-5,1-k*1e-5] for k in range(11)])
res = solve_ivp(multi_Lorenz,[1,25],u0, method="LSODA")

plt.figure(figsize=(14,6))
plt.plot(res.t, res.y[::3].T)
plt.grid(); plt.show()

The visible split of the solution family is at t=22, however in the differences to the k=0 graph as reference, the increase to the tenfold of the initial difference happens almost immediately in the spike around t=1 and remains in that range up to t=12

N = max(k for k,t in enumerate(res.t) if t<12)
plt.figure(figsize=(14,6))
plt.plot(res.t[:N], (res.y[3::3,:N]-res.y[0,:N]).T)
plt.grid(); plt.show()