I am currently trying to solve a second order BVP, but I do not think I am handling breaking down the second order differential equation into two first order equations while taking into account the complex term correctly.

The differential equation is
f”(z) = -2ikf'(z) + V(phi_sol(z))f(z)

Boundary Conditions
z-> inf,
f'(z) = 0, f(z) = 1

I am trying to determine the small fluctuations about a soliton solution that I have determined earlier and phi_sol(z) is that solution of the soliton at z and k is the wavenumber/mode of the fluctuations.

Code (in order to run the code you may have to pip install pynverse):

import numpy as np
from pynverse import inversefunc
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt

def f(x, a=1):
    term1 = (x**2 - 1) * (x**2 + a**2)
    term2 = -2 * np.arctan(x / a) + a * (np.log(1 - x) - np.log(1 + x))
    numerator = term1 * term2
    denominator = 2 * np.sqrt(2) * (a + a**2) * np.sqrt((x**2 - 1)**2 * 
                  (x**2 + a**2)**2)
    return numerator / denominator

# Set up the domain
domain = np.linspace(-0.999, 0.999, 500)

# Evaluate the function to determine its range (image)
image_values = f(domain)


# Define the inverse of the function using pynverse within the valid 
# image range
inverse_f = inversefunc(f, domain=[-0.999, 0.999])

# Restrict y_values to the valid range of the function
y_values = np.linspace(np.min(image_values), np.max(image_values), 500)


a = 1.0  
z = y_values


# Define the second derivative of the potential U''(phi)
def U_prime_prime(phi):
   return 4 * (14 * phi**6 + (15 * a**2 - 15) * phi**4 + (3 * a**4 - 
              12 * a**2 + 3) * phi**2 - a**4 + a**2)

# Define the effective potential V_eff(z)
def V_eff(z,  U_prime_prime):
    # We need to compute V_eff based on the soliton profile 
    # phi_sol(z)
    return U_prime_prime(inverse_f(z)) - U_prime_prime(1)
    # U_prime_prime(1) is a constant

def system(z, y, k,  U_prime_prime):
    # y[0] = f(z), y[1] = f'(z)
    f, fp = y
    Veff = V_eff(z,  U_prime_prime)
    # First order system of equations
    return np.vstack([fp, -2j * k * fp + Veff * f])

# Boundary conditions
def bc(ya, yb):
    # Boundary conditions: f(z) -> 1 as z -> ∞, f'(z) -> 0 as z -> ∞
    return np.array([ya[0], yb[1]-1])

# Initial guess for f(z)
y_init = np.zeros((2, z.size))
y_init[0] = np.cos(z)  

k = 1.0  

solution = solve_bvp(lambda z, y: system(z, y, k,  U_prime_prime), bc, z, y_init)

if solution.success:
    # Plot the solution f(z)
    plt.plot(z, solution.sol(z)[0], label="f(z)")
    plt.xlabel("z")
    plt.ylabel("f(z)")
    plt.title("Solution to the Quantum Fluctuation Equation")
    plt.legend()
    plt.grid(True)
    plt.show()
else:
    print("Solution failed to converge")