def METHOD_NAME(dim, n, order):
    # ### problem definition
    import dolfin as df
    if dim == 2:
        mesh = df.UnitSquareMesh(n, n)
    elif dim == 3:
        mesh = df.UnitCubeMesh(n, n, n)
    else:
        raise NotImplementedError
    V = df.FunctionSpace(mesh, 'CG', order)
    g = df.Constant(1.0)
    c = df.Constant(1.)
    class DirichletBoundary(df.SubDomain):
        def inside(self, x, on_boundary):
            return abs(x[0] - 1.0) < df.DOLFIN_EPS and on_boundary
    db = DirichletBoundary()
    bc = df.DirichletBC(V, g, db)
    u = df.Function(V)
    v = df.TestFunction(V)
    f = df.Expression('x[0]*sin(x[1])', degree=2)
    F = df.inner((1 + c*u**2)*df.grad(u), df.grad(v))*df.dx - f*v*df.dx
    df.solve(F == 0, u, bc,
             solver_parameters={'newton_solver': {'relative_tolerance': 1e-6}})
    # ### pyMOR wrapping
    from pymor.bindings.fenics import FenicsOperator, FenicsVectorSpace, FenicsVisualizer
    from pymor.models.basic import StationaryModel
    from pymor.operators.constructions import VectorOperator
    space = FenicsVectorSpace(V)
    op = FenicsOperator(F, space, space, u, (bc,),
                        parameter_setter=lambda mu: c.assign(mu['c'].item()),
                        parameters={'c': 1},
                        solver_options={'inverse': {'type': 'newton', 'rtol': 1e-6}})
    rhs = VectorOperator(op.range.zeros())
    fom = StationaryModel(op, rhs,
                          visualizer=FenicsVisualizer(space))
    return fom