Python and model solving

I am currently try to solve a relatively simple model. The core equations attached in the following file:

Equations

I am writing V_U and V_B only as functions of d for convenience, but they obviously depend on many other parameters. When trying to solve for it, my code looks as follows:

def V_G(SMIC,beta,sigma):
return u(SMIC,sigma)/(1-beta)

def V_U(b,sigma,beta,RSA,w_B,SMIC,p,d,dmax):
if d==0:
    return u(RSA,sigma)
else:
    return u(b,sigma) + beta*(p*max(V_G(SMIC,beta, 
    sigma),V_U(b,sigma,beta,RSA,w_B,SMIC,p,d- 
    1,dmax),V_B(b,sigma,beta,w_B,SMIC,p,d-1,dmax)) + (1- 
    p)*max(V_U(b,sigma,beta,RSA,w_B,SMIC,p,d- 
    1,dmax),V_B(b,sigma,beta,w_B,SMIC,p,d-1,dmax)))

def V_B(b,sigma,beta,w_B,SMIC,p,d,dmax):
if d==dmax:
    return u(w_B,sigma) +  beta*max(V_B(b,sigma,beta,w_B,SMIC,p,d,dmax),V_U(b,sigma,beta,RSA,w_B,SMIC,p,d,dmax))
else:
    return u(w_B,sigma) + beta*max(V_B(b,sigma,beta,w_B,SMIC,p,d+1,dmax),V_U(b,sigma,beta,RSA,w_B,SMIC,p,d,dmax))


Unemp= V_U(650,0.2,0.997,500,w_B,1500,0.1,1,24)    

I unfortunately get stuck in some sort of infinite loop. I am not familiar with python, and my math skills are quite limited. I don't know whether this model cannot actually be numerically simulated, or if the method I am using is just not appropriate. Maybe I should use an equation solver ?

Any help, comment or advice would be much appreciated.