def func(t,y,a,A,w):  # тут твоя функция идет
    return np.array( [a*(-((y[0]**3)/3 - y[0]) + y[1]),
                      -1*y[0] + A*cos(w*t)])

def nordsieck_func(z0,z1,t,h,a,A,w): #эта хуита идет в fsolve

    M = np.array([[1, 1-2/3, 1-4/3],
                  [0,     0,     0],
                  [0,  -1/3, 1-2/3]])
    l = np.array([[2/3,1,1/3]]).transpose()
    z0 = z0.reshape((3,2))
    z1 = z1.reshape((3,2))

    return (h*l*np.array([[a*(-((z0[0][0]**3)/3 - z0[0][0]) + z0[0][1]), -1*z0[0][0] + A*cos(w*t)]]) + np.dot(M,z1) - z0).flatten()

def nordsieck(a,A,w,T, d = (0.01,0.001), falloff = 0.01, graphs = {1,2,3,4,5}, dim = 5, ignore_hcalc = False):
    y = np.array([[2,0]])
    t = 0
    h = 0.001
    t_arr = [0]
    h_arr = []
    z = np.array([y[-1], h*func(t,y[-1],a,A,w), np.array([0,0])])

    calculation = True
    while calculation:

        dmax_buffer = 0
        h_calculation = True
        while h_calculation:

            flat_z = z.flatten()
            z_new_flat = fsolve(nordsieck_func,flat_z,(flat_z,t,h,a,A,w))
            z_new = z_new_flat.reshape((3,2))
            if ignore_hcalc:
                break

            delta = z_new[0] - z[0]
            d_max = np.amax(np.absolute(delta))

            if abs(d_max - dmax_buffer) > falloff: # на быстрых участках нужно экспоненциально менять шаг для
                dmax_buffer = d_max                # линейного уменьшения прироста
            else:                                  # если находим такой участок, перестаем менять шаг
                h_calculation = False

            if d_max > d[0]:
                h -= h*0.1
            elif d_max < d[1]:
                h += h*0.1
            else:
                h_calculation = False

        z = z_new
        t += h
        h_arr.append(h)
        y = np.append(y, [z_new[0]], axis = 0)
        t_arr.append(t)

        if t > T:
            calculation = False

    if graphs:
        n = 1
        m1 = ceil(np.sqrt(len(graphs)))
        m2 = ceil(len(graphs)/m1)
        fig = plt.figure(figsize = (m2*dim,m1*dim))

        if 1 in graphs:
            ax = fig.add_subplot(m1, m2, n)
            ax.plot(t_arr[1:],h_arr)
            n+=1

        if 2 in graphs:
            ax = fig.add_subplot(m1, m2, n)
            ax.plot(y.transpose()[0],y.transpose()[1])
            n+=1
        if 3 in graphs:
            ax = fig.add_subplot(m1, m2, n)
            ax.plot(t_arr,y.transpose()[0])
            n+=1
        if 4 in graphs:
            ax = fig.add_subplot(m1, m2, n)
            ax.plot(t_arr,y.transpose()[1])
            n+=1

        if 5 in graphs:
            ax = fig.add_subplot(m1, m2, n, projection='3d')
            ax.plot(y.transpose()[0],y.transpose()[1],t_arr)
            n+=1

        fig.show()
    return [y,t_arr]