Python: how to interpret y from solve_ivp

Question

This is surely a trivial question, but it prevents my complete understanding of solve_ivp from scypy.integrate on which I am presently training ... The image below defines the problem I'm trying to solve with solve_ivp:

So, in order to find y(t), I specify the function to integrate, the initial values, the time span, and then I run solve_ivp, as shown in the code below:

# Function to integrate
def fun(t, u):    
    x1 = u[0]     # "u": function to found / 4 components x1, x2, x3 and x4
    x2 = u[1]
    x3 = u[2]
    x4 = u[3]

    dx1_dt = 1   # "u'(t) = F(t,u(t))": derivatives of components
    dx2_dt = x3
    dx3_dt = x4
    dx4_dt = np.exp(x1) + 5*x2 - x1*x3
    
    return [dx1_dt, dx2_dt, dx3_dt, dx4_dt]

# Specify initial conditions
x1_0 = 0.0
x2_0 = 0.0
x3_0 = 0.0
x4_0 = 0.0

y_0 = np.array([x1_0, x2_0, x3_0, x4_0])

# Specify initial and final times
t0 = 0.0  
tf = 10.0 

t_span = np.array([t0, tf])

# Resolution
position = solve_ivp(fun, t_span, y_0, method='RK45', max_step=0.1)

Now, solve_ivp returns a ndarray named y (in the present case, it will be position.y, of shape (4, 104)) which according to the scipy.integrate.solve_ivp documentation, gives the "Values of the solution at t".

So far, it's good.

My Question is just:

In the present problem, what are the values giving y(t) : y[0], y[1], y[2] or y[3]? As far as I understood how solve_ivp works, it should be y[1], corresponding to the second line of vector u(t). Is it right?

Answer 1

Yes, that is correct, position.y[1] contains the values for the solution function y(t). You should also find that position.y[0] coincides with position.t.

If you want additional values or a faster computation, use the teval or dense_output options. With the max_step argument you are forcing more than 100 internal steps. Without that and one of the other options, the number of internal steps is adapted to the internal error tolerance and usually less. You get the values at the desired points then from the piecewise polynomial interpolation function, either implicitly with teval=... or explicitly with the inclusion of the "dense output" interpolation function in the return object as position.sol.