Runge-Kutta method overflow encountered in double_scalars

Question

I'm no programming expert, but I'm trying to model a runge-kutta method basically by trial an error of a boundary condition. And that's fine by me except that if a try to expand the time domain (i.e add more items to the vector by increasing xStop), the program suddenly explodes and overflows.

Can you tell my why is it so unstable?

import numpy as np
from matplotlib import pyplot as plt

def runge_kutta4(F,x,y,xStop,h):

    def rk4(F,x,y,h):
        K0 = h*F(x,y)
        K1 = h*F(x + h/2.0, y + K0/2.0)
        K2 = h*F(x + h/2.0, y + K1/2.0)
        K3 = h*F(x + h, y + K2)
        return (K0 + 2.0*K1 + 2.0*K2 + K3)/6.0
    X = []
    Y = []
    X.append(x)
    Y.append(y)
    while x < xStop:
        h = min(h,xStop - x)
        y = y + rk4(F,x,y,h)
        x = x + h
        X.append(x)
        Y.append(y)
    return np.array(X),np.array(Y)

def secant(f,x1,x2,tol,mx):
    for i in range(mx):
        xnew = x2-(x2-x1)/(f(x2)-f(x1)) * f(x2)
        if abs(xnew-x2) < tol: break
        else:
            x1 = x2
            x2 = xnew
    else:
        print('cagasteZ')
    return xnew


def initCond(u): # Initial values of [y,y’,y"];
# use 'u' if unknown
    return np.array([0.0, 0.0, u])

beta = 7.5 #CONDICION DE BORDE DEL OTRO EXTREMO

def res(u): # Boundary condition residual--see Eq. (8.3)
    X,Y = runge_kutta4(F,xStart,initCond(u),xStop,h)
    y = Y[len(Y) - 1]
    residuo = y[0] - beta   #ESTE ES EL VALOR DEL OTRO EXTREMO
    return residuo

def F(x,y): # First-order differential equations
    F = np.zeros(3)
    F[0] = y[1]
    F[1] = y[2]
    F[2] = y[1]**2-y[0]*y[2]-1 
    return F

xStart = 0.0 # Start of integration
xStop = 7.5 # End of integration
u1 = 1 # 1st trial value of unknown init. cond.
u2 = 1.5 # 2nd trial value of unknown init. cond.
h = 0.01 # initial step size
u = secant(res,u1,u2,0.001,1000)#
X,Y = runge_kutta4(F,xStart,initCond(u),xStop,h)
plt.figure()
plt.plot(X,Y[:,0],'-',X,Y[:,1],'-',X,Y[:,2],'-')
plt.xlabel('$\eta$')
plt.legend(('f','df/d$\eta$','d$^2$f/d$\eta^2$'),loc = 2)
plt.grid(True)
plt.show()