<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def getPoly(rCut, nMax):
 """Used to calculate discrete vectors for the polynomial basis functions. Args: rCut(float):
 Radial cutoff nMax(int):
 Number of polynomial radial functions """

    rCutVeryHard = rCut+5.0
    rx = 0.5*rCutVeryHard*(x + 1)

    basisFunctions = []
    for i in range(1, nMax + 1):
        basisFunctions.append(lambda rr, i=i, rCut=rCut: (rCut - np.clip(rr, 0, rCut))**(i+2))

    # Calculate the overlap of the different polynomial functions in a
    # matrix S. These overlaps defined through the dot product over the
    # radial coordinate are analytically calculable: Integrate[(rc - r)^(a
    # + 2) (rc - r)^(b + 2) r^2, {r, 0, rc}]. Then the weights B that make
    # the basis orthonormal are given by B=S^{-1/2}
    S = np.zeros((nMax, nMax))
    for i in range(1, nMax+1):
        for j in range(1, nMax+1):
            S[i-1, j-1] = (2*(rCut)**(7+i+j))/((5+i+j)*(6+i+j)*(7+i+j))
    betas = sqrtm(np.linalg.inv(S))

    # If the result is complex, the calculation is currently halted.
    if (betas.dtype == np.complex128):
        raise ValueError(
            "Could not calculate normalization factors for the polynomial basis"
            " in the domain of real numbers. Lowering the number of radial "
            "basis functions is advised."
        )

    fs = np.zeros([nMax, len(x)])
    for n in range(1, nMax+1):
        fs[n-1, :] = (rCut-np.clip(rx, 0, rCut))**(n+2)

    gss = np.dot(betas, fs)

    return nMax, rx, gss