How to use the dymos.RungeKutta function in dymos

def make_brachistochrone_phase(transcription='gauss-lobatto', num_segments=8, transcription_order=3,
                               compressed=True):

    if transcription == 'gauss-lobatto':
        t = dm.GaussLobatto(num_segments=num_segments,
                            order=transcription_order,
                            compressed=compressed)
    elif transcription == 'radau-ps':
        t = dm.Radau(num_segments=num_segments,
                     order=transcription_order,
                     compressed=compressed)
    elif transcription == 'runge-kutta':
        t = dm.RungeKutta(num_segments=num_segments,
                          order=transcription_order,
                          compressed=compressed)

    phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)

    return phase

#
        # Initialize the Problem and the optimization driver
        #
        p = om.Problem(model=om.Group())
        p.driver = om.ScipyOptimizeDriver()
        p.driver.declare_coloring()

        #
        # Create a trajectory and add a phase to it
        #
        traj = p.model.add_subsystem('traj', dm.Trajectory())

        phase = traj.add_phase('phase0',
                               dm.Phase(ode_class=BrachistochroneODE,
                                        transcription=dm.RungeKutta(num_segments=10)))

        #
        # Set the variables
        #
        phase.set_time_options(initial_bounds=(0, 0), duration_bounds=(.5, 10))

        phase.add_state('x', rate_source=BrachistochroneODE.states['x']['rate_source'],
                        units=BrachistochroneODE.states['x']['units'],
                        fix_initial=True, fix_final=False, solve_segments=False)

        phase.add_state('y', rate_source=BrachistochroneODE.states['y']['rate_source'],
                        units=BrachistochroneODE.states['y']['units'],
                        fix_initial=True, fix_final=False, solve_segments=False)

        phase.add_state('v', rate_source=BrachistochroneODE.states['v']['rate_source'],
                        targets=BrachistochroneODE.states['v']['targets'],

# define a Trajectory object and add to model
    traj = dm.Trajectory()
    p.model.add_subsystem('traj', subsys=traj)

    # define a Transcription
    if transcription == 'gauss-lobatto':
        t = dm.GaussLobatto(num_segments=num_segments,
                            order=transcription_order,
                            compressed=compressed)
    elif transcription == 'radau-ps':
        t = dm.Radau(num_segments=num_segments,
                     order=transcription_order,
                     compressed=compressed)
    elif transcription == 'runge-kutta':
        t = dm.RungeKutta(num_segments=num_segments,
                          order=transcription_order,
                          compressed=compressed)

    # define a Phase as specified above and add to Phase
    if not delay:
        phase = dm.Phase(ode_class=vanderpol_ode, transcription=t)
    else:
        phase = dm.Phase(ode_class=vanderpol_ode_group, transcription=t)  # distributed component group
    traj.add_phase(name='phase0', phase=phase)

    t_final = 15.0
    phase.set_time_options(fix_initial=True, fix_duration=True, duration_val=t_final, units='s')

    # set the State time options
    phase.add_state('x0', fix_initial=False, fix_final=False,
                    rate_source='x0dot',

p.driver.options['optimizer'] = optimizer
    p.driver.options['dynamic_simul_derivs'] = True

    if optimizer == 'SNOPT':
        p.driver.opt_settings['Major iterations limit'] = 1000
        p.driver.opt_settings['iSumm'] = 6
        p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-6
        p.driver.opt_settings['Major optimality tolerance'] = 1.0E-6
        p.driver.opt_settings['Function precision'] = 1.0E-12
        p.driver.opt_settings['Linesearch tolerance'] = 0.1
        p.driver.opt_settings['Major step limit'] = 0.5
        # p.driver.opt_settings['Verify level'] = 3

    t = {'gauss-lobatto': GaussLobatto(num_segments=num_seg, order=transcription_order),
         'radau-ps': Radau(num_segments=num_seg, order=transcription_order),
         'runge-kutta': RungeKutta(num_segments=num_seg)}

    phase = Phase(ode_class=MinTimeClimbODE, transcription=t[transcription])

    p.model.add_subsystem('phase0', phase)

    phase.set_time_options(fix_initial=True, duration_bounds=(50, 400),
                           duration_ref=100.0)

    phase.set_state_options('r', fix_initial=True, lower=0, upper=1.0E6,
                            ref=1.0E3, defect_ref=1.0E3, units='m')

    phase.set_state_options('h', fix_initial=True, lower=0, upper=20000.0,
                            ref=1.0E2, defect_ref=1.0E2, units='m')

    phase.set_state_options('v', fix_initial=True, lower=10.0,
                            ref=1.0E2, defect_ref=1.0E2, units='m/s')

p = Problem(model=Group())

    if optimizer == 'SNOPT':
        p.driver = pyOptSparseDriver()
        p.driver.options['optimizer'] = optimizer
        p.driver.opt_settings['Major iterations limit'] = 100
        p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-6
        p.driver.opt_settings['Major optimality tolerance'] = 1.0E-6
        p.driver.opt_settings['iSumm'] = 6
    else:
        p.driver = ScipyOptimizeDriver()

    p.driver.options['dynamic_simul_derivs'] = dynamic_simul_derivs

    if transcription == 'runge-kutta':
        transcription = RungeKutta(num_segments=num_segments, compressed=compressed)
    elif transcription == 'gauss-lobatto':
        transcription = GaussLobatto(num_segments=num_segments,
                                     order=transcription_order,
                                     compressed=compressed)
    elif transcription == 'radau-ps':
        transcription = Radau(num_segments=num_segments,
                              order=transcription_order,
                              compressed=compressed)

    phase = Phase(ode_class=BrachistochroneVectorStatesODE,
                  transcription=transcription)

    p.model.add_subsystem('phase0', phase)

    phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))

OPT, OPTIMIZER = set_pyoptsparse_opt(optimizer, fallback=True)
    p.driver.options['optimizer'] = OPTIMIZER

    if simul_derivs:
        p.driver.declare_coloring()

    if transcription == 'gauss-lobatto':
        t = dm.GaussLobatto(num_segments=num_segments,
                            order=transcription_order,
                            compressed=compressed)
    elif transcription == 'radau-ps':
        t = dm.Radau(num_segments=num_segments,
                     order=transcription_order,
                     compressed=compressed)
    elif transcription == 'runge-kutta':
        t = dm.RungeKutta(num_segments=num_segments,
                          order=transcription_order,
                          compressed=compressed)

    phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)

    p.model.add_subsystem('phase0', phase)

    phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10), units='s')

    phase.add_state('x', fix_initial=True, fix_final=False, solve_segments=False,
                    units='m', rate_source='xdot')
    phase.add_state('y', fix_initial=True, fix_final=False, solve_segments=False,
                    units='m', rate_source='ydot')
    phase.add_state('v', fix_initial=True, fix_final=False, solve_segments=False,
                    units='m/s', rate_source='vdot', targets=['v'])

def make_traj(transcription='gauss-lobatto', transcription_order=3, compressed=True,
              connected=False):

    t = {'gauss-lobatto': GaussLobatto(num_segments=20, order=transcription_order, compressed=compressed),
         'radau': Radau(num_segments=20, order=transcription_order, compressed=compressed),
         'runge-kutta': RungeKutta(num_segments=20, compressed=compressed)}

    traj = Trajectory()

    traj.add_design_parameter('c', opt=False, val=1.5, units='DU/TU')

    # First Phase (burn)

    burn1 = Phase(ode_class=FiniteBurnODE, transcription=t[transcription])

    burn1 = traj.add_phase('burn1', burn1)

    burn1.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
    burn1.set_state_options('r', fix_initial=True, fix_final=False, defect_scaler=100.0)
    burn1.set_state_options('theta', fix_initial=True, fix_final=False, defect_scaler=100.0)
    burn1.set_state_options('vr', fix_initial=True, fix_final=False, defect_scaler=100.0)
    burn1.set_state_options('vt', fix_initial=True, fix_final=False, defect_scaler=100.0)