def METHOD_NAME(a0_a_ratio, sh, sv, nu, a0, G, initialFrictionAngle, finalFrictionAngle, param_m):
	a = a0/a0_a_ratio
	xi_well = 1.0 - a0/a # the auxiliary variable xi at the wellbore, xi = (a-a0)/a = 1-1/(a/a0)
	nPoints = 1000
	initialFrictionAngle *= np.pi/180.0 # converted to rad
	finalFrictionAngle *= np.pi/180.0 # converted to rad
	param_b_i = compute_param_b(initialFrictionAngle)
	param_b_f = compute_param_b(finalFrictionAngle)
	# ELastic trial
	pw = 2.0*G*xi_well + sh
	r_e,sr_e,s0_e,sz_e = epwAnal.solution_elastic(sh,sv,1.0,pw)
	p_e, q_e = epwAnal.invariants(sr_e[0],s0_e[0],sz_e[0])
	if ((q_e/p_e)<param_b_i): # Pure elastic wellbore
		return [0],[0],[0],[0],r_e,sr_e,s0_e,sz_e,pw,p_e,q_e
	else: # Elastic-plastic wellbore
		q_ep_boundary = compute_q_ep_boudary_EDP(sh,sv,param_b_i)
		# Elastic-Plastic boundary
		sr_ep_boundary, s0_ep_boundary, sz_ep_boundary = elastic_plastic_boundary_EDP(sh,sv,q_ep_boundary)
		sr_p = [sr_ep_boundary]
		s0_p = [s0_ep_boundary]
		sz_p = [sz_ep_boundary]
		epsV_p = [0] # strain calculated starting from the reference state that is the ep boundary
		# Eq. 36 in Chen and Abousleiman 2017
		xi_ep_boundary = (sr_ep_boundary - sh) / 2.0 / G
		# Mesh from elastic-plastic boundary to wellbore
		dxi = (xi_well - xi_ep_boundary) / (nPoints - 1)
		xi = np.linspace(xi_ep_boundary, xi_well, nPoints)
		for i in range(1, nPoints):
			p,q = epwAnal.invariants(sr_p[i-1],s0_p[i-1],sz_p[i-1])
			exp_epsilon_V = np.exp(epsV_p[i-1])# epsV is the volumetric strain from the elastic-plastic boundary
			E = 2.0 * G * (1.0 + nu)
			param_b = q/p # yield condition F=0
			dFdp, dFdq, dFdbeta = dFdpq_EDP(p,q)
			dBetadGammaP = hardening_EDP(param_b,param_m,param_b_f,param_b_i)
			dGdq = dFdq # assicated model
			dFdHardningParam,dHardningParamdPlasticVar,dPotentialdStressVar = dFdbeta,dBetadGammaP,dGdq
			sr_i, s0_i, sz_i, epsV_i = epwAnal.solution_plastic(sr_p[i - 1], s0_p[i - 1], sz_p[i - 1], epsV_p[i-1], xi[i - 1], dxi, dFdp, dFdq, E, nu, dFdHardningParam,dHardningParamdPlasticVar,dPotentialdStressVar)
			sr_p.append(sr_i)
			s0_p.append(s0_i)
			sz_p.append(sz_i)        
			epsV_p.append(epsV_i)
		# Wellbore surface stress
		pw = sr_i
		p = (sr_i + s0_i + sz_i)/3.0
		q = np.sqrt(0.5) * np.sqrt( (sr_i-s0_i)**2.0 + (sr_i-sz_i)**2.0 + (s0_i-sz_i)**2.0 )
		# Compute the normalized radial coordinate r/a
		r_p = epwAnal.compute_radialCoordinate(xi, epsV_p)
		# Elastic zone
		r_ep_boundary = r_p[0]
		r_e,sr_e,s0_e,sz_e = epwAnal.solution_elastic(sh,sv,r_ep_boundary,sr_ep_boundary)
		return r_p,sr_p,s0_p,sz_p,r_e,sr_e,s0_e,sz_e,pw,p,q