How to use the dymos.utils.lgl.lgl function in dymos

import openmdao.api as om
        from openmdao.utils.assert_utils import assert_near_equal

        import dymos as dm
        from dymos.examples.battery_multibranch.battery_multibranch_ode import BatteryODE
        from dymos.utils.lgl import lgl

        prob = om.Problem()

        opt = prob.driver = om.ScipyOptimizeDriver()
        opt.declare_coloring()
        opt.options['optimizer'] = 'SLSQP'

        num_seg = 5
        seg_ends, _ = lgl(num_seg + 1)

        traj = prob.model.add_subsystem('traj', dm.Trajectory())

        # First phase: normal operation.
        transcription = dm.Radau(num_segments=num_seg, order=5, segment_ends=seg_ends, compressed=False)
        phase0 = dm.Phase(ode_class=BatteryODE, transcription=transcription)
        traj_p0 = traj.add_phase('phase0', phase0)

        traj_p0.set_time_options(fix_initial=True, fix_duration=True)
        traj_p0.add_state('state_of_charge', fix_initial=True, fix_final=False,
                          targets=['SOC'], rate_source='dXdt:SOC')

        # Second phase: normal operation.

        phase1 = dm.Phase(ode_class=BatteryODE, transcription=transcription)
        traj_p1 = traj.add_phase('phase1', phase1)

import openmdao.api as om
        from openmdao.utils.assert_utils import assert_near_equal

        import dymos as dm

        from dymos.examples.aircraft_steady_flight.aircraft_ode import AircraftODE
        from dymos.examples.plotting import plot_results
        from dymos.utils.lgl import lgl

        p = om.Problem(model=om.Group())
        p.driver = om.pyOptSparseDriver()
        p.driver.options['optimizer'] = 'SLSQP'
        p.driver.declare_coloring()

        num_seg = 15
        seg_ends, _ = lgl(num_seg + 1)

        traj = p.model.add_subsystem('traj', dm.Trajectory())

        phase = traj.add_phase('phase0',
                               dm.Phase(ode_class=AircraftODE,
                                        transcription=dm.Radau(num_segments=num_seg,
                                                               segment_ends=seg_ends,
                                                               order=3, compressed=False)))

        # Pass Reference Area from an external source
        assumptions = p.model.add_subsystem('assumptions', om.IndepVarComp())
        assumptions.add_output('S', val=427.8, units='m**2')
        assumptions.add_output('mass_empty', val=1.0, units='kg')
        assumptions.add_output('mass_payload', val=1.0, units='kg')

        phase.set_time_options(initial_bounds=(0, 0),

p = Problem(model=Group())
    p.driver = pyOptSparseDriver()
    _, optimizer = set_pyoptsparse_opt(optimizer, fallback=False)
    p.driver.options['optimizer'] = optimizer
    p.driver.options['dynamic_simul_derivs'] = True
    if optimizer == 'SNOPT':
        p.driver.opt_settings['Major iterations limit'] = 20
        p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-6
        p.driver.opt_settings['Major optimality tolerance'] = 1.0E-6
        p.driver.opt_settings["Linesearch tolerance"] = 0.10
        p.driver.opt_settings['iSumm'] = 6
    if optimizer == 'SLSQP':
        p.driver.opt_settings['MAXIT'] = 50

    num_seg = 15
    seg_ends, _ = lgl(num_seg + 1)

    phase = Phase(ode_class=AircraftODE,
                  transcription=Radau(num_segments=num_seg, segment_ends=seg_ends,
                                      order=3, compressed=compressed, solve_segments=solve_segments))

    # Pass Reference Area from an external source
    assumptions = p.model.add_subsystem('assumptions', IndepVarComp())
    assumptions.add_output('S', val=427.8, units='m**2')
    assumptions.add_output('mass_empty', val=1.0, units='kg')
    assumptions.add_output('mass_payload', val=1.0, units='kg')

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

    phase.set_time_options(initial_bounds=(0, 0),
                           duration_bounds=(300, 10000),
                           duration_ref=5600)

self.initial_state_vec = np.zeros(self.state_vec_size)

        # Setup the control interpolants
        if self.options['control_options']:
            for name, options in self.options['control_options'].items():
                self.add_input(name='controls:{0}'.format(name),
                               val=np.ones(((ncdsps,) + options['shape'])),
                               units=options['units'],
                               desc='Values of control {0} at control discretization '
                                    'nodes within the segment.'.format(name))
                interp = LagrangeBarycentricInterpolant(control_disc_seg_stau, options['shape'])
                self.options['ode_integration_interface'].control_interpolants[name] = interp

        if self.options['polynomial_control_options']:
            for name, options in self.options['polynomial_control_options'].items():
                poly_control_disc_ptau, _ = lgl(options['order'] + 1)
                self.add_input(name='polynomial_controls:{0}'.format(name),
                               val=np.ones(((options['order'] + 1,) + options['shape'])),
                               units=options['units'],
                               desc='Values of polynomial control {0} at control discretization '
                                    'nodes within the phase.'.format(name))
                interp = LagrangeBarycentricInterpolant(poly_control_disc_ptau, options['shape'])
                self.options['ode_integration_interface'].polynomial_control_interpolants[name] = \
                    interp

        if self.options['design_parameter_options']:
            for name, options in self.options['design_parameter_options'].items():
                self.add_input(name='design_parameters:{0}'.format(name), val=np.ones(options['shape']),
                               units=options['units'],
                               desc='values of design parameter {0}.'.format(name))

        if self.options['input_parameter_options']:

self.rate_jacs = {}
        self.rate2_jacs = {}
        self.val_jac_rows = {}
        self.val_jac_cols = {}
        self.rate_jac_rows = {}
        self.rate_jac_cols = {}
        self.rate2_jac_rows = {}
        self.rate2_jac_cols = {}
        self.sizes = {}

        self.add_input('t_duration', val=1.0, units=self.options['time_units'],
                       desc='duration of the phase to which this interpolated control group '
                            'belongs')

        for name, options in self.options['polynomial_control_options'].items():
            disc_nodes, _ = lgl(options['order'] + 1)
            num_control_input_nodes = len(disc_nodes)
            shape = options['shape']
            size = np.prod(shape)
            units = options['units']
            rate_units = get_rate_units(units, self.options['time_units'], deriv=1)
            rate2_units = get_rate_units(units, self.options['time_units'], deriv=2)

            input_shape = (num_control_input_nodes,) + shape
            output_shape = (num_nodes,) + shape

            L_de, D_de = lagrange_matrices(disc_nodes, eval_nodes)
            _, D_dd = lagrange_matrices(disc_nodes, disc_nodes)
            D2_de = np.dot(D_de, D_dd)

            self._matrices[name] = L_de, D_de, D2_de

subsets = {
        'disc': np.arange(0, n, 2, dtype=int),
        'state_disc': np.zeros((1,), dtype=int) if seg_idx == 0 or shooting == 'multiple'
        else np.empty((0,), dtype=int),
        'state_input': np.zeros((1,), dtype=int) if seg_idx == 0 or shooting == 'multiple'
        else np.empty((0,), dtype=int),
        'control_disc': np.arange(n, dtype=int),
        'control_input': np.arange(n, dtype=int) if not compressed or seg_idx == 0
        else np.arange(1, n, dtype=int),
        'segment_ends': np.array([0, n-1], dtype=int),
        'col': np.arange(1, n, 2, dtype=int),
        'all': np.arange(n, dtype=int),
        'solution': np.arange(n, dtype=int),
    }

    return subsets, lgl(n)[0]

else np.arange(1, n, dtype=int),
        'segment_ends': np.array([0, n-1], dtype=int),
        'col': np.arange(1, n, 2, dtype=int),
        'all': np.arange(n, dtype=int),
        'solution': np.arange(n, dtype=int),
    }

    subsets['solver_solved'] = subsets['state_input'] if compressed or \
        seg_idx == 0 else subsets['state_input'][1::]

    idxs_not_in_solved = np.where(np.in1d(subsets['state_input'],
                                          subsets['solver_solved'],
                                          invert=True))[0]
    subsets['solver_indep'] = subsets['state_input'][idxs_not_in_solved]

    return subsets, lgl(n)[0]

self.rate_jacs = {}
        self.rate2_jacs = {}
        self.val_jac_rows = {}
        self.val_jac_cols = {}
        self.rate_jac_rows = {}
        self.rate_jac_cols = {}
        self.rate2_jac_rows = {}
        self.rate2_jac_cols = {}
        self.sizes = {}

        self.add_input('t_duration', val=1.0, units=self.options['time_units'],
                       desc='duration of the phase to which this interpolated control group '
                            'belongs')

        for name, options in self.options['polynomial_control_options'].items():
            disc_nodes, _ = lgl(options['order'] + 1)
            num_control_input_nodes = len(disc_nodes)
            shape = options['shape']
            size = np.prod(shape)
            units = options['units']
            rate_units = get_rate_units(units, self.options['time_units'], deriv=1)
            rate2_units = get_rate_units(units, self.options['time_units'], deriv=2)

            input_shape = (num_control_input_nodes,) + shape
            output_shape = (num_output_nodes,) + shape

            L_do, D_do = lagrange_matrices(disc_nodes, output_nodes_ptau)
            _, D_dd = lagrange_matrices(disc_nodes, disc_nodes)
            D2_do = np.dot(D_do, D_dd)

            self._matrices[name] = L_do, D_do, D2_do