def f(x, y):
    return x + y


def eiler(f, a, b, y0, n):
    h = (b - a) / n
    x = [k for k in range(n + 1)]
    y = [k for k in range(n + 1)]
    x[0] = a
    y[0] = y0
    for i in range(n):
        x[i + 1] = x[i] + h
        ypr = y[i] + h / 2 * f(x[i], y[i])
        y[i + 1] = y[i] + h * f(x[i] + h / 2, ypr)

    return x, y


def kut(f, a, b, y0, n):
    h = (b - a) / n
    x = [k for k in range(n + 1)]
    y = [k for k in range(n + 1)]
    x[0] = a
    y[0] = y0
    for i in range(n):
        x[i + 1] = x[i] + h
        k1 = h * f(x[i], y[i])
        k2 = h * f(x[i] + h / 2, y[i] + k1 / 2)
        k3 = h * f(x[i] + h / 2, y[i] + k2 / 2)
        k4 = h * f(x[i] + h, y[i] + k3)
        y[i + 1] = y[i] + 1 / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
    return x, y


def met_adam_bosh(fun, a, b, y0, n):
    h = (b - a) / n
    x = [x for x in range(n + 1)]
    y = [x for x in range(n + 1)]
    f = [x for x in range(n + 1)]
    x[0] = a
    y[0] = y0
    f[0] = fun(x[0], y[0])
    for i in range(3):
        x[i + 1] = x[i] + h
        k1 = h * f[i]
        k2 = h * fun(x[i] + h / 2, y[i] + k1 / 2)
        k3 = h * fun(x[i] + h / 2, y[i] + k2 / 2)
        k4 = h * fun(x[i] + h, y[i] + k3)
        y[i + 1] = y[i] + 1 / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
        f[i + 1] = fun(x[i + 1], y[i + 1])
    for i in range(4, n):
        x[i + 1] = x[i] + h
        y[i + 1] = y[i] + h / 24 * (55 * f[i] - 59 * f[i - 1] +
                                    37 * f[i - 2] - 9 * f[i - 3])
        f[i + 1] = fun(x[i + 1], y[i + 1])
    return x, y


def met_adam_multon(fun, a, b, y0, n):
    h = (b - a) / n
    x = [x for x in range(n + 1)]
    y = [x for x in range(n + 1)]
    f = [x for x in range(n + 1)]
    x[0] = a
    y[0] = y0
    f[0] = fun(x[0], y[0])
    for i in range(3):
        x[i + 1] = x[i] + h
        k1 = h * f[i]
        k2 = h * fun(x[i] + h / 2, y[i] + k1 / 2)
        k3 = h * fun(x[i] + h / 2, y[i] + k2 / 2)
        k4 = h * fun(x[i] + h, y[i] + k3)
        y[i + 1] = y[i] + 1 / 6 * (k1 + 2 * k2 + 2 * k3 + k4)
        f[i + 1] = fun(x[i + 1], y[i + 1])
    for i in range(4, n):
        x[i + 1] = x[i] + h
        ypred = y[i] + h / 24 * (55 * f[i] - 59 * f[i - 1] +
                                 37 * f[i - 2] - 9 * f[i - 3])
        fpred = fun(x[i + 1], ypred)
        y[i + 1] = y[i] + h / 24 * (9 * fpred + 19 * f[i] - 5 * f[i - 1] + f[i - 2])
        f[i + 1] = fun(x[i + 1], y[i + 1])
    return x, y


def main():
    a = 0
    b = 5
    y0 = 1
    n = 20
    x, y = eiler(f, a, b, y0, n)
    print(x)
    print(y)

    x, y = kut(f, a, b, y0, n)
    print(x)
    print(y)

    x, y = met_adam_multon(f, a, b, y0, n)
    print(x)
    print(y)

    x, y = met_adam_bosh(f, a, b, y0, n)
    print(x)
    print(y)