How to get the derivative of a variable with respect to time at a point of the solution of a differential equation solved with solve_ivp

I wrote some code solving a simulation (with scipy.solve_ivp, lsoda method) of a physical system. Inside this simulation I am running a control loop to control a physical aspect of the system. This control loop uses a variable that calculates the error between a reference system and the current state of the physical system. Problem for the control loop is that I need to get the integral and derivative of that error to calculate the output of the control loop. While I can simply get the integral of this error by including that in my dy vector of the solve_ivp function (and getting the y output) I have no idea how to get the derivative. Is there a way to do this while still using the lsoda method of the solve_ivp function?

I have added a code snippet to show an abstraction of the problem. It would be even nicer if this could be solved by simulating something akin to; derivative: newest_var - old_var / dt (with set intervals of constant dt) inside the differential solver since our control system in real life will use this and having a way to measure the error induced by using this relatively innacurate method on the physical system will ensure we have a good idea of what our total error could be irl.

import numpy as np

def rhs_func(t, rhs_integrated):
    pos = rhs_integrated[0]
    vel = rhs_integrated[1]
    var_error_integrated = rhs_integrated[2]

    # these are the gains for the control, in this case they are completely random
    i_gain = 0.1  # integral gain
    p_gain = 0.1  # proportional gain
    d_gain = 0.1  # derivative gain
    # suppose you know what you want the pos to be at a certain velocity
    pos_control = np.exp(0.25*vel)
    # with this you can compute the error between what you want it to be 
    # and what it needs to be
    var_error = pos - pos_control
    # now here lies the problem I have a derivative gain but no idea how to get the 
    # derivative of var error (the part I commented out)
    vel = rhs_integrated[1] - p_gain*var_error + i_gain*var_error_integrated # + d_gain*var_error_derivative
    acc = 10 - 0.5*vel  
    return np.hstack((vel, acc, var_error))  # dy vector

# suppose the starting velocity is 10 acc is 0 and var error is also 0
sol = solve_ivp(rhs_func, [0, 10], [10, 0, 0], t_eval=[1,2,3,4,5,6,7,8,9,10])
print(sol)