題名の通り、\vec{z}'=fという2階の微分方程式(数値的に解くために一階のベクトルの微分方程式の形にしている)を、初期値X0とパラメータepsilonを変化させて、無限遠で0になるという境界条件に合う(x0,epsilon)の組み合わせを見つけて、その時のepsilonの最小値を拾ってくるということをpythonでしたいのですが、以下のようなコードを書いたところ、やはりループがいけないのかとっても時間がかかってしまいます。本当はmeshgridなどで出来れば一番早いと思うのですが、うまくいきませんでした。

極力ループを使わないような方法で以下のコードと同じようなことができる方法があればご教授お願い致します。

def sol(q_x,kappa,h,x0,epsilon):

    y = Symbol('y')

    def fsol(x,y,epsilon):
        return f(x,y,q_x,kappa,h,epsilon)

    y = np.linspace(0.0,1000,10001)

    return odeint(fsol,x0,y,args=(epsilon,))

alpha = np.linspace(0,2*np.pi,101)
X = np.cos(alpha)
Y = np.sin(alpha)
ZEROS = [[0.0]*2 for i in range(len(alpha))]
X0 = np.c_[X,Y,ZEROS]

def dispersion0(q_x,kappa,h,cutoff):

    def z(x0,epsilon):
        return sol(q_x,kappa,h,x0,epsilon)

    epsilon = np.linspace(0,q_x**2+kappa+h,101)

    Phi1 = np.empty([alpha.shape[0],epsilon.shape[0]])
    Phi2 = np.empty([alpha.shape[0],epsilon.shape[0]])
    for i in range(0,len(alpha)):
        for j in range(0,len(epsilon)):
            Phi1[i][j] = z(X0[i],epsilon[j])[1000][0]
            Phi2[i][j] = z(X0[i],epsilon[j])[1000][1]
    energy = 0
    for k in range(0,len(epsilon)):
        if abs(Phi1[:,k])<cutoff & abs(Phi2[:,k])<cutoff :
            energy = epsilon(k)
            break            

    return energy