def METHOD_NAME():
    # dims
    num_mass = 4
    nx = 2*num_mass
    nu = num_mass-1
    # symbolic variables
    sym_x = ca.SX.sym('x', nx, 1) # states
    sym_u = ca.SX.sym('u', nu, 1) # controls
    sym_xdot = ca.SX.sym('xdot', casadi_length(sym_x)) #state derivatives
    # dynamics
    # continuous time
    Ac = np.zeros((nx, nx))
    for ii in range(num_mass):
        Ac[ii,num_mass+ii] = 1.0
        Ac[num_mass+ii,ii] = -2.0
    for ii in range(num_mass-1):
        Ac[num_mass+ii,ii+1] = 1.0
        Ac[num_mass+ii+1,ii] = 1.0
    Bc = np.zeros((nx, nu))
    for ii in range(nu):
        Bc[num_mass+ii, ii] = 1.0
    c_const = np.zeros(nx)
    # discrete time
    Ts = 0.5 # sampling time
    M = scipy.linalg.expm(np.vstack(( np.hstack((Ts*Ac, Ts*Bc)), np.zeros((nu, int(2*nx/2+nu))))))
    A = M[:nx,:nx]
    B = M[:nx,nx:]
    expr_f_expl = Ac@sym_x + Bc@sym_u + c_const
    expr_f_impl = expr_f_expl - sym_xdot
    expr_phi = A@sym_x + B@sym_u
    # constraints
    expr_h = ca.vertcat(sym_u, sym_x)
    expr_h_e = sym_x
    # nonlnear least squares
    expr_y = ca.vertcat(sym_u, sym_x)
    expr_y_e = sym_x
    # external cost
    yr_u = np.zeros(nu)
    yr_x = np.zeros(nx)
    dWu = 2*np.ones(nu)
    dWx = np.ones(nx)
    ymyr = ca.vertcat(sym_u, sym_x) - ca.vertcat(yr_u, yr_x)
    ymyr_e = sym_x - yr_x
    expr_ext_cost = 0.5 * (ymyr.T @ ( np.concatenate((dWu, dWx)) *  ymyr))
    expr_ext_cost_e = 0.5 * ymyr_e.T @ (dWx * ymyr_e)
    # populate structure
    model = {}
    model['nx'] = nx
    model['nu'] = nu
    model['sym_x'] = sym_x
    model['sym_xdot'] = sym_xdot
    model['sym_u'] = sym_u
    model['expr_f_expl'] = expr_f_expl
    model['expr_f_impl'] = expr_f_impl
    model['expr_phi'] = expr_phi
    model['expr_h'] = expr_h
    model['expr_h_e'] = expr_h_e
    model['expr_y'] = expr_y
    model['expr_y_e'] = expr_y_e
    model['expr_ext_cost'] = expr_ext_cost
    model['expr_ext_cost_e'] = expr_ext_cost_e
    return model