ctrl, clsys = ct.create_statefbk_iosystem(
      sys, ([g1, ..., gN], [p1, ..., pN]), gainsched_indices=[s1, ..., sq])

    # Create a simple nonlinear system to check (kinematic car)
    def unicycle_update(t, x, u, params):
        return np.array([np.cos(x[2]) * u[0], np.sin(x[2]) * u[0], u[1]])

    def unicycle_output(t, x, u, params):
        return x

    unicycle = ct.NonlinearIOSystem(
        unicycle_update, unicycle_output, inputs = ['v', 'phi'], 
        outputs = ['x', 'y', 'theta'], states = ['x_', 'y_', 'theta_'])

    # Speeds and angles at which to compute the gains
    speeds = [1, 5, 10]
    angles = np.linspace(0, pi/2, 4)
    points = list(itertools.product(speeds, angles))

    # Gains for each speed (using LQR controller)
    Q = np.identity(unicycle.nstates)
    R = np.identity(unicycle.ninputs)
    gains = [np.array(ct.lqr(unicycle.linearize(
        [0, 0, angle], [speed, 0]), Q, R)[0]) for speed, angle in points]

    #
    # Schedule on desired speed and angle
    #

    # Create gain scheduled controller
    ctrl, clsys = ct.create_statefbk_iosystem(
        unicycle, (gains, points), gainsched_indices=['vd', 'theta'])