Problem with recursive implementation of Subspace Iteration method in Numpy

Question

I am having trouble with implementing the method of subspace iteration to find the eigenvalues and vectors of a random, symmetric matrix, A that is mxm with m = 10. The function that I have written to perform this task is

def subspaceiterate(A,V,v,j):
    if j == 0:
        return v
    else:
        v_jm1 = V[:,j-1]
        v_jm1 = np.reshape(v_jm1,(np.size(V,axis=0),1))
        v = v - np.matmul(v_jm1.T,np.matmul(A,v_jm1))
        j = j - 1
        return subspaceiterate(A,V,v,j)

where A is my symmetric data matrix, V is an mxm matrix whose columns represent the current, calculated state of the eigenvectors, v is a particular mx1 column ofV that represents a particular eigenvector and j is the position of vin V.

As I understand it, the basic idea of Subspace Iteration is to simultaneously perform a Power Iteration on the different v's while subtracting away the components of v that are in the direction of the eigenvectors associated with eigenvalues larger than the one associated with the particular v in question. I have tried to design this function to operate in such a capacity but the results I get do not support the aforementioned idea. All the eigenvalues converge quickly (as they should), but do so to a value very close to the largest one instead of to the different, distinct eigenvalues (whose values I know because I have calculated them beforehand with np.linalg.eigs()).

The v's are normalized in the program before and after the function is called and the eigenvalues calculated with the Rayleigh-Taylor quotient v_j.T*A*v_j (implemented with np.matmul()). At this point I am lost and hoping someone brilliant can shine some light and help me out. `