def METHOD_NAME(self, c=0):
    r"""
    Return a vector `x` such that ``self(x) == c``.
    INPUT:
    - ``c`` -- (default: 0) a rational number
    OUTPUT: A non-zero vector `x` satisfying ``self(x) == c``.
    ALGORITHM:
    Uses PARI's :pari:`qfsolve`. Algorithm described by Jeroen Demeyer; see comments on :trac:`19112`
    EXAMPLES::
        sage: F = DiagonalQuadraticForm(QQ, [1, -1]); F
        Quadratic form in 2 variables over Rational Field with coefficients:
        [ 1 0 ]
        [ * -1 ]
        sage: F.solve()
        (1, 1)
        sage: F.solve(1)
        (1, 0)
        sage: F.solve(2)
        (3/2, -1/2)
        sage: F.solve(3)
        (2, -1)
    ::
        sage: F = DiagonalQuadraticForm(QQ, [1, 1, 1, 1])
        sage: F.solve(7)
        (1, 2, -1, -1)
        sage: F.solve()
        Traceback (most recent call last):
        ...
        ArithmeticError: no solution found (local obstruction at -1)
    ::
        sage: Q = QuadraticForm(QQ, 2, [17, 94, 130])
        sage: x = Q.solve(5); x
        (17, -6)
        sage: Q(x)
        5
        sage: Q.solve(6)
        Traceback (most recent call last):
        ...
        ArithmeticError: no solution found (local obstruction at 3)
        sage: G = DiagonalQuadraticForm(QQ, [5, -3, -2])
        sage: x = G.solve(10); x
        (3/2, -1/2, 1/2)
        sage: G(x)
        10
        sage: F = DiagonalQuadraticForm(QQ, [1, -4])
        sage: x = F.solve(); x
        (2, 1)
        sage: F(x)
        0
    ::
        sage: F = QuadraticForm(QQ, 4, [0, 0, 1, 0, 0, 0, 1, 0, 0, 0]); F
        Quadratic form in 4 variables over Rational Field with coefficients:
        [ 0 0 1 0 ]
        [ * 0 0 1 ]
        [ * * 0 0 ]
        [ * * * 0 ]
        sage: F.solve(23)
        (23, 0, 1, 0)
    Other fields besides the rationals are currently not supported::
        sage: F = DiagonalQuadraticForm(GF(11), [1, 1])
        sage: F.solve()
        Traceback (most recent call last):
        ...
        TypeError: solving quadratic forms is only implemented over QQ
    """
    if self.base_ring() is not QQ:
        raise TypeError("solving quadratic forms is only implemented over QQ")
    M = self.Gram_matrix()
    # If no argument passed for c, we just pass self into qfsolve().
    if not c:
        x = qfsolve(M)
        if isinstance(x, Integer):
            raise ArithmeticError("no solution found (local obstruction at {})".format(x))
        return x
    # If c != 0, define a new quadratic form Q = self - c*z^2
    d = self.dim()
    N = Matrix(self.base_ring(), d+1, d+1)
    for i in range(d):
        for j in range(d):
            N[i, j] = M[i, j]
    N[d, d] = -c
    # Find a solution x to Q(x) = 0, using qfsolve()
    x = qfsolve(N)
    # Raise an error if qfsolve() doesn't find a solution
    if isinstance(x, Integer):
        raise ArithmeticError("no solution found (local obstruction at {})".format(x))
    # Let z be the last term of x, and remove z from x
    z = x[-1]
    x = x[:-1]
    # If z != 0, then Q(x/z) = c
    if z:
        return x * (1/z)
    # Case 2: We found a solution self(x) = 0. Let e be any vector such
    # that B(x,e) != 0, where B is the bilinear form corresponding to self.
    # To find e, just try all unit vectors (0,..0,1,0...0).
    # Let a = (c - self(e))/(2B(x,e)) and let y = e + a*x.
    # Then self(y) = B(e + a*x, e + a*x) = self(e) + 2B(e, a*x)
    #              = self(e) + 2([c - self(e)]/[2B(x,e)]) * B(x,e) = c.
    e = vector([1] + [0] * (d-1))
    i = 0
    while self.bilinear_map(x, e) == 0:
        e[i] = 0
        i += 1
        e[i] = 1
    a = (c - self(e)) / (2 * self.bilinear_map(x, e))
    return e + a*x