<p dir="auto">I am trying to solve a simple example with the <code class="notranslate">dopri5</code> integrator in <code class="notranslate">scipy.integrate.ode</code>. As the documentation states</p>
<blockquote>
<p dir="auto">This is an explicit runge-kutta method of order (4)5 due to Dormand &amp; Prince (with stepsize control and dense output).</p>
</blockquote>
<p dir="auto">this should work. So here is my example:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt


def MassSpring_with_force(t, state, f):
    &quot;&quot;&quot; Simple 1DOF dynamics model: m ddx(t) + k x(t) = f(t)&quot;&quot;&quot;
    # unpack the state vector
    x = state[0]
    xd = state[1]

    # these are our constants
    k = 2.5 # Newtons per metre
    m = 1.5 # Kilograms

    # compute acceleration xdd
    xdd = ( ( -k*x + f) / m )

    # return the two state derivatives
    return [xd, xdd]


def force(t):
    &quot;&quot;&quot; Excitation force &quot;&quot;&quot;
    f0 = 1  # force amplitude [N]
    freq = 20  # frequency[Hz]
    omega = 2 * np.pi *freq  # angular frequency [rad/s]
    return f0 * np.sin(omega*t)


# Time range
t_start = 0
t_final = 1

# Main program
state_ode_f = ode(MassSpring_with_force)

state_ode_f.set_integrator('dopri5', rtol=1e-4, nsteps=500, 
                           first_step=1e-6, max_step=1e-1, verbosity=True)

state2 = [0.0, 0.0]  # initial conditions
state_ode_f.set_initial_value(state2, 0)
state_ode_f.set_f_params(force(0))

sol = np.array([[t_start, state2[0], state2[1]]], dtype=float)

print(&quot;Time\t\t Timestep\t dx\t\t ddx\t\t state_ode_f.successful()&quot;)

while state_ode_f.successful() and state_ode_f.t &lt; (t_final):
    state_ode_f.set_f_params(force(state_ode_f.t))
    state_ode_f.integrate(t_final, step=True)
    sol = np.append(sol, [[state_ode_f.t, state_ode_f.y[0], state_ode_f.y[1]]], axis=0)
    print(&quot;{0:0.8f}\t {1:0.4e} \t{2:10.3e}\t {3:0.3e}\t {4}&quot;.format(
            state_ode_f.t, sol[-1, 0]- sol[-2, 0], state_ode_f.y[0], state_ode_f.y[1], state_ode_f.successful()))"><pre class="notranslate"><code class="notranslate">import numpy as np
from scipy.integrate import ode
import matplotlib.pyplot as plt


def MassSpring_with_force(t, state, f):
    """ Simple 1DOF dynamics model: m ddx(t) + k x(t) = f(t)"""
    # unpack the state vector
    x = state[0]
    xd = state[1]

    # these are our constants
    k = 2.5 # Newtons per metre
    m = 1.5 # Kilograms

    # compute acceleration xdd
    xdd = ( ( -k*x + f) / m )

    # return the two state derivatives
    return [xd, xdd]


def force(t):
    """ Excitation force """
    f0 = 1  # force amplitude [N]
    freq = 20  # frequency[Hz]
    omega = 2 * np.pi *freq  # angular frequency [rad/s]
    return f0 * np.sin(omega*t)


# Time range
t_start = 0
t_final = 1

# Main program
state_ode_f = ode(MassSpring_with_force)

state_ode_f.set_integrator('dopri5', rtol=1e-4, nsteps=500, 
                           first_step=1e-6, max_step=1e-1, verbosity=True)

state2 = [0.0, 0.0]  # initial conditions
state_ode_f.set_initial_value(state2, 0)
state_ode_f.set_f_params(force(0))

sol = np.array([[t_start, state2[0], state2[1]]], dtype=float)

print("Time\t\t Timestep\t dx\t\t ddx\t\t state_ode_f.successful()")

while state_ode_f.successful() and state_ode_f.t &lt; (t_final):
    state_ode_f.set_f_params(force(state_ode_f.t))
    state_ode_f.integrate(t_final, step=True)
    sol = np.append(sol, [[state_ode_f.t, state_ode_f.y[0], state_ode_f.y[1]]], axis=0)
    print("{0:0.8f}\t {1:0.4e} \t{2:10.3e}\t {3:0.3e}\t {4}".format(
            state_ode_f.t, sol[-1, 0]- sol[-2, 0], state_ode_f.y[0], state_ode_f.y[1], state_ode_f.successful()))
</code></pre></div>
<p dir="auto">The result I get is:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="Time         Timestep        dx          ddx         state_ode_f.successful()
1.00000000   1.0000e+00      0.000e+00   0.000e+00   True"><pre class="notranslate"><code class="notranslate">Time         Timestep        dx          ddx         state_ode_f.successful()
1.00000000   1.0000e+00      0.000e+00   0.000e+00   True
</code></pre></div>
<p dir="auto">Hence, only one time-step is computed which is obviously incorrect.</p>
<p dir="auto">This works with <code class="notranslate">vode</code> and <code class="notranslate">zvode</code> integrators</p>
<p dir="auto"><em>Original ticket <a href="http://projects.scipy.org/scipy/ticket/1820" rel="nofollow">http://projects.scipy.org/scipy/ticket/1820</a> on 2013-01-22 by trac user timdiller, assigned to trac user teoliphant.</em></p>
<p dir="auto">The dopri5 and dop853 methods of scipy.integrate.ode do not support the step keyword argument. This is nicely illustrated in this [http://stackoverflow.com/questions/12926393/using-adaptive-step-sizes-with-scipy-integrate-ode StackOverflow article].</p>
<p dir="auto">Implementing this would provide a better-working alternative to MATLAB's ode45(), etc.</p>