def Ufunc(n, Q, V, P):
    """
    The U operator, see Messiah book, chapter 16.

    Parameters:
    ----------
    n: int
        Order of expansion.
    Q: matrix
        Projection matrix rescaled by inverse energy difference.
    V: matrix
        Perturbing part of Hamiltonian.
    P: matrix
        Projection matrix to degenerate subspace we are perturbing from.

    Returns:
    --------
    matrix
        The expansion of U to nth order.
    """
    if n == 0:
        return P
    elif n > 0:
        U = V*Ufunc(n-1, Q, V, P)
        for p in range(1,n):
            U = U - Ufunc(p,Q,V,P)*V*Ufunc(n-p-1,Q,V,P)
        return Q*U