import numpy as np

def rk2(t, y, h, f):
    k_1 = h*f(t,y)
    k_2 = h*f(t + h/2.0, y+k_1/2.0)
    y = y + k_2
    return y

def rk4(t, y, h, f):
    k_1 = h*f(t,y)
    k_2 = h*f(t + h/2.0, y + k_1/2.0)
    k_3 = h*f(t + h/2.0, y + k_2/2.0)
    k_4 = h*f(t + h, y + k_3)
    y = y + (1.0/6.0)*(k_1 + 2.0*k_2 + 2.0*k_3 +k_4)
    return y

def rk45(t, y, h, f, err_tol, h_min, h_max):
    k_1 = h*f(t, y)
    k_2 = h*f(t + h/4.0, y + k_1/4.0)
    k_3 = h*f(t + 3.0*h/8.0, y + 3.0*k_1/32.0 +9.0*k_2/32.0)
    k_4 = h*f(t + 12.0*h/13.0, y+1932.0*k_1/2197.0 - 7200.0*k_2/2197.0 +7296.0*k_3/2197.0)
    k_5 = h*f(t + h, y + 439.0*k_1/216.0 - 8.0*k_2 + 3680.0*k_3/513.0 -845.0*k_4/4104.0)
    k_6 = h*f(t + h/2.0, y - 8.0*k_1/27.0 + 2.0*k_2 - 3544.0*k_3/2565.0 + 1859.0*k_4/4104.0 - 11.0*k_5/40.0) 
    Y = y + 25.0*k_1/216.0 + 1408.0*k_3/2565.0 + 2197.0*k_4/4101.0 - k_5/5.0
    Z = y + 16.0*k_1/135.0 + 6656.0*k_3/12825.0 +28561.0*k_4/56430.0 - 9.0*k_5/50.0 + 2.0*k_6/55.0
    err = np.abs(Z - Y)
    if (err[0] == 0 or err[1] == 0):
        s = 0.0
    else:
        s = 0.84*np.power(err_tol*h/err[0], 0.25)
    if (s < 0.75 and h > 2*h_min):
        h /= 2.0
    elif (s > 1.5 and h < h_max/2.0):
        h *= 2.0
    return h,Z
        
