def RK4_singlestep(b, theta0tpm, omega0tpm, t, dt, rhs, theta0, q, omegad):
    """ take a single RK-4 timestep from t to t+dt for the system 
        ydot = rhs """

    theta = theta0tpm
    omega = omega0tpm
    
    # get the RHS at several points
    thetadot1, omegadot1 = rhs(b, theta, omega,t, q, omegad)
    thetadot2, omegadot2 = rhs(b, theta+0.5*dt*thetadot1, omega+0.5*dt*omegadot1,t, q, omegad)
    thetadot3, omegadot3 = rhs(b, theta+0.5*dt*thetadot2, omega+0.5*dt*omegadot2,t, q, omegad)
    thetadot4, omegadot4 = rhs(b, theta+0.5*dt*thetadot3, omega+0.5*dt*omegadot3,t, q, omegad)
            
    # advance
    thetanew = theta + (dt/6.0)*(thetadot1 + 2.0*thetadot2 + 2.0*thetadot3 + thetadot4)
    omeganew = omega + (dt/6.0)*(omegadot1 + 2.0*omegadot2 + 2.0*omegadot3 + omegadot4)

    return thetanew, omeganew