I am trying to find the solution of DEQ dy/dx = 10*e^(-88.8888889(x-2)^2) - 0.6y.

Function:

def dydx(x,y):
    return -0.6*y + 10*np.exp(-(x-2)**2/(2*.075**2))

Single Step Solution:

def direct(x,y,h):

    k1 = dydx(x,y)
    k2 = dydx(x+h/2, y+k1*h/2)
    k3 = dydx(x+h/2, y+k2*h/2)
    k4 = dydx(x+h, y+k3*h)

    y = y + (1/6)*(k1+2*k2+2*k3+k4)*h

    return y

Double Step Solution:

def halving(x,y,h):

    h = h/2
    xh = np.zeros(3)
    yh = np.zeros(3)

    xh[0] = x
    yh[0] = y

    for j in range(0,2):
    
        k1 = dydx(xh[j],yh[j])
        k2 = dydx(xh[j]+h/2, yh[j]+k1*h/2)
        k3 = dydx(xh[j]+h/2, yh[j]+k2*h/2)
        k4 = dydx(xh[j]+h, yh[j]+k3*h)

        yh[j+1] = yh[j] + (1/6)*(k1+2*k2+2*k3+k4)*h
        xh[j+1] = xh[j] + h

    return yh[2]

Find the adaptive step size 'h'

def adapt(x,y,h):
    error = abs(direct(x,y,h)) - abs(halving(x,y,h))  # it checks whether h falls between a range or not

    if error <= 0.01 and error >= .001:
        return h
    elif error > 0.01:
        return adapt(x,y,h/2) # if error is large recalling adapt function with h/2
    elif error < .001:
        return adapt(x,y,2*h) # 2*h if error is small 

Main Driver Program

def rk4(xi,yi,h,xlim):

    nval = int((xlim-xi)/h)

    x = np.zeros(nval+1)
    y = np.zeros(nval+1)

    x[0] = xi
    y[0] = yi

    for i in range(0,nval):
    
        h = adapt(x[i],y[i],h)
        y[i+1] = direct(x[i],y[i],h)
        x[i+1] = x[i] + h 

    return x,y

Evaluate using

rk4(0,0.5,.1,4)

I am using jupyter notebook. This program crashes. Probably the recursion section. But why? The error should have been within range after a few recursive call?