def METHOD_NAME():
    tape = Tape()
    set_working_tape(tape)
    mesh = UnitSquareMesh(10, 10)
    V_element = VectorElement("CG", mesh.ufl_cell(), 1)
    V = FunctionSpace(mesh, V_element)
    u = TrialFunction(V)
    v = TestFunction(V)
    W = FunctionSpace(mesh, "CG", 2)
    c = project(Expression("sin(x[0])*sin(x[1])", degree=1, annotate=False), W, annotate=False)
    expr = Expression(("c*c*c", "2*c*c"), c=c, degree=2)
    first_deriv = Expression(("3*c*c", "4*c"), c=c, degree=2, annotate=False)
    second_deriv = Expression(("6*c", "4"), c=c, degree=2, annotate=False)
    expr.user_defined_derivatives = {c: first_deriv}
    first_deriv.user_defined_derivatives = {c: second_deriv}
    a = inner(u,v)*dx
    L = expr**2*inner(expr, v)*dx
    bc = DirichletBC(V, (0, 0), "on_boundary")
    u_ = Function(V)
    solve(a == L, u_, bc)
    J = assemble(inner(u_, u_)*dx)
    control = Control(c)
    Jhat = ReducedFunctional(J, control)
    h = Function(c.function_space())
    h.vector()[:] = rand(c.function_space().dim())*4.0
    dJdm = h._ad_dot(Jhat.derivative())
    Hm = h._ad_dot(Jhat.hessian(h))
    assert(taylor_test(Jhat, c, h, dJdm=dJdm, Hm=Hm) > 2.9)