def vcv_recursive(putative_root, vcv_matrix, finished):
    """
    This is a fast recursive function to calculate the variance co-variance matrix. 
    
    Input(s):
    putative_root - the root node (initially) and subsequently, other clades from a Bio.Phylo tree object
    vcv_matrix - the matrix which will store everything. Needs to be n x n where n is the 
            total number of terminals in the tree object where putative_root belongs. Should be zeroes to start
    finished - a list of all called nodes to keep track of everything for record keeping

    Output(s):
    vcv_matrix - the updated vcv_matrix described above
    finished - the updated finished list described above
    """
    terminals = putative_root.get_terminals()
    if not set(terminals).issubset(set(finished)):
        vcv_matrix[len(finished):len(finished)+len(terminals), len(finished):len(finished)+len(terminals)] += putative_root.branch_length
    if len(putative_root.clades) == 2:
            vcv_matrix, finished = vcv_recursive(putative_root.clades[0], vcv_matrix, finished)
            vcv_matrix, finished = vcv_recursive(putative_root.clades[1], vcv_matrix, finished)
    elif len(putative_root.clades) == 0:
        finished.append(putative_root)
    return vcv_matrix, finished