def METHOD_NAME(T=50):
    n_x = 500
    nt = int(T / 0.01) + 1
    wave_speed = 0.1
    l = 1.
    dx = l / (n_x-1)
    # construct H_op
    space = NumpyVectorSpace(n_x)
    Dxx = diags(
        [-2 * np.ones(n_x), np.ones(n_x-1), np.ones(n_x-1)],
        [0, 1, -1],
        format='csr',
    )
    H_op = BlockDiagonalOperator([
        NumpyMatrixOperator(-wave_speed**2 / dx * Dxx),
        LincombOperator([IdentityOperator(space)], [1/dx]),
    ])
    # construct initial_data
    h = lambda s: (0 <= s) * (s <= 1) * (1 - 3/2 * s**2 + 3/4 * s**3) + (1 < s) * (s <= 2) * ((2-s)**3)/4
    bump = lambda xi: h(np.abs(4*(xi - l/2)))
    initial_data = H_op.source.make_array([
        space.make_array(bump(np.linspace(0, l, n_x))),
        space.make_array(np.zeros(n_x)),
    ])
    fom = QuadraticHamiltonianModel(T, initial_data, H_op, nt=nt, name='hamiltonian_wave_equation')
    # fom.operator = fom.operator.with_(solver_options={'type': 'to_matrix'}) #TODO
    return fom