W-cycle multigrid Python

Question

I'm having troubles on solving a PDE by means of multigrid method in Python. I particular I have to implement a W-cycle.

My function takes in input $uh$ as a starting vector, defined as follows, and the mesh size is 1024.

The inputs are:

def w_i(i, x): return np.sin(i * x * np.pi)

n= 1024

x = np.linspace(0, 1, n+1)

uh = 0.1 * w_i(1, x) + 0.5 * w_i(n//2, x) + w_i(n-3, x)

fh = np.zeros_like(x)

What i have done is a recoursive formula:

def w_cycle(uh,fh):

N = len(uh) - 1
h = 1/N
if N == 2:
    cvec = c(2)
    uh[1] = fh[1]/((2/h**2)+c[1]) 
    smax = 0
else:
    for i in range(alpha1):
        jacobi(fh, uh)    #pre-smoothing
    Ah = CreateA(N)  #it creates the iteration matrix on the grid corrsponding to h
    f2h = I_h_2h(fh - (Ah @ uh))   #I_h_2h is a restriction operator
    u2h = np.zeros_like(f2h)
    u2h = w_cycle(u2h,f2h)
    u2h = w_cycle(u2h,f2h)
    uh = uh + I_2h_h(u2h)          #I_2h_h is a prolongation operator
    for i in range(alpha2):
        save = uh.copy()
        jacobi(fh, uh)             #post-smoothing
    smax = np.linalg.norm(uh-save,np.inf)
return uh

I should obtain by executing 10 times this program that the norm infinity of the resulting $uh$ must be of order of 10^-10 but what i obtain is of order of 10^-2 so there must be some errors in the code. In particular I think that $uh$ is not always updated correctly, while I'm sure that all the other functions that are in the code are correct. Can someone help me to underestand where I'm wrong?