Untitled

from scipy import optimize as opt 
import numpy as np
from scipy.optimize import minimize
import matplotlib.pyplot as plt 
from scipy.integrate import quad
from scipy.misc import derivative
cb = 223.6
cf = 14.8
pi = np.pi
m_h = 125
m_w = 80.385
m_z = 91.1876
m_t = 173.07
v = 246.22
M = 600
mw2 = lambda phi: (m_w**2)*(phi**2/v**2)
mz2 = lambda phi: (m_z**2)*(phi**2/v**2)
mt2 = lambda phi: (m_t**2)*(phi**2/v**2)
    
def jb(y_t):
  return quad((lambda x: x**2 * (np.log(1-np.exp(-np.sqrt(x**2+y_t**2))))), 0, np.infty)[0]
def jf(z_t):
  return quad((lambda x: x**2 * (np.log(1+np.exp(-np.sqrt(x**2+z_t**2))))), 0, np.infty)[0]
 
   
def V(phi,T, M):
    A = 3*v**4/(4*M**2) - m_h**2/2
    B = -3*v**2/(2*M**2) + m_h**2/(2*v**2)
    V0 = A*(phi**2)/2 + (B*phi**4)/4 + (1/(8*M**2))*phi**6  
    MASSA_W = ((mw2(phi)**2)*np.log(mw2(phi)/m_w**2)) -((3*mw2(phi)**2)/2) + (2*mw2(phi))*(m_w**2) - (m_w**4)/2
    MASSA_Z = ((mz2(phi)**2)*np.log(mz2(phi)/m_z**2)) -((3*mz2(phi)**2)/2) + (2*mz2(phi))*(m_z**2) - (m_z**4)/2
    MASSA_t = ((mt2(phi)**2)*np.log(mt2(phi)/m_t**2)) -((3*mt2(phi)**2)/2) + (2*mt2(phi))*(m_t**2) - (m_t**4)/2
    V1 = (1/(64*(pi**2)))*(6*MASSA_W+3*MASSA_Z-12*MASSA_t)
    if T!=0:
        del_V1 = ((T**4/(2*pi**2))*((6*jb(np.sqrt(mw2(phi))/T) + (3*jb(np.sqrt(mz2(phi))/T)) - (12*jf(np.sqrt(mt2(phi))/T)))))         
        return V0+V1+del_V1
    else:
        return V0+V1
def v_linha(phi,T):
    return V(phi,T, M)-V(1e-15,T, M)


y0 = 0
rf = 10
d_V = derivative(V,1.2,dx=1e-6)


def f(r,s):
    x,w = s
    dx_dr = w
    dv_dr = -2*w/r - d_V
    return dx_dr, dv_dr

def yf(w0):
    sol = integrate.solve_ivp(f,(0,rf),(y0,w0))
    x,w = sol.x
    return x[-1]

yf_vec = np.vectorize(yf)
w0 = optimize.newton(yf_vec,10)
r = np.linspace(0,rf,100)
sol = integrate.solve_ivp(f, (0,rf),(y0,w0),t_eval=r)
x,v = sol.y
plt.plot(r,x,".")
print(w0)