def RK4(f, y_i, t_i, dt):
    """Performs Runge-Kutta integration.

    Parameters
    ----------
    f : callable
        Derivative function.

    y_i : ndarray
        Current state.

    t_i : float
        Current time.

    dt : float
        Time step.

    Returns
    -------
    y_j : ndarray
        State at next time step.
    """

    # Determine k constants
    k0 = f(y_i, t_i)
    k1 = f(y_i+0.5*dt*k0, t_i+0.5*dt)
    k2 = f(y_i+0.5*dt*k1, t_i+0.5*dt)
    k3 = f(y_i+dt*k2, t_i+dt)

    # Calculate y at next timestep
    return y_i+0.16666666666666666666667*(k0+2.0*k1+2.0*k2+k3)*dt