Changes

hessian_matrix = array([[0, -2*sin(theta)*cos(theta), 0, 0, 0, 0], [-2*sin(theta)*cos(theta), 2*R1*sin(theta)**2 - 2*R1*cos(theta)**2 - 2*(R1rho_p + kex*phi_ex/(kex**2 + we**2))*sin(theta)**2 + 2*(R1rho_p + kex*phi_ex/(kex**2 + we**2))*cos(theta)**2, 2*sin(theta)*cos(theta), 2*kex*sin(theta)*cos(theta)/(kex**2 + we**2), 2*(-2*kex**2*phi_ex/(kex**2 + we**2)**2 + phi_ex/(kex**2 + we**2))*sin(theta)*cos(theta), -4*kex*phi_ex*we*sin(theta)*cos(theta)/(kex**2 + we**2)**2], [0, 2*sin(theta)*cos(theta), 0, 0, 0, 0], [0, 2*kex*sin(theta)*cos(theta)/(kex**2 + we**2), 0, 0, -2*kex**2*sin(theta)**2/(kex**2 + we**2)**2 + sin(theta)**2/(kex**2 + we**2), -2*kex*we*sin(theta)**2/(kex**2 + we**2)**2], [0, 2*(-2*kex**2*phi_ex/(kex**2 + we**2)**2 + phi_ex/(kex**2 + we**2))*sin(theta)*cos(theta), 0, -2*kex**2*sin(theta)**2/(kex**2 + we**2)**2 + sin(theta)**2/(kex**2 + we**2), (8*kex**3*phi_ex/(kex**2 + we**2)**3 - 6*kex*phi_ex/(kex**2 + we**2)**2)*sin(theta)**2, (8*kex**2*phi_ex*we/(kex**2 + we**2)**3 - 2*phi_ex*we/(kex**2 + we**2)**2)*sin(theta)**2], [0, -4*kex*phi_ex*we*sin(theta)*cos(theta)/(kex**2 + we**2)**2, 0, -2*kex*we*sin(theta)**2/(kex**2 + we**2)**2, (8*kex**2*phi_ex*we/(kex**2 + we**2)**3 - 2*phi_ex*we/(kex**2 + we**2)**2)*sin(theta)**2, 8*kex*phi_ex*we**2*sin(theta)**2/(kex**2 + we**2)**3 - 2*kex*phi_ex*sin(theta)**2/(kex**2 + we**2)**2]])