def f(xi, yi):
    ys = xi + yi
    return ys


# Gegeben:
x0 = 0
y0 = 0
h = 0.2
xmax = 0.4
Ey = 0.0001

x = [x0]
y = [[y0]]

n = 0
while x[-1] < xmax:
    x.append(x[-1]+h)
    fn = f(x[n], y[n][-1])
    if n == 0:
        y_neu = y[n][-1] + h * fn
    else:
        y_neu = y[n-1][-1] + 2*h * fn
    y.append([y_neu])
    dy = 1
    while dy > Ey:
        y_neu = y[n][-1] + h/2 * (fn + f(x[n+1], y[n+1][-1]))
        y[-1].append(y_neu)
        dy = abs(y[-1][-2]-y[-1][-1])
    n += 1

print(x)
print(y)

werte = []
for i in range(len(x)):
    werte.append((x[i], y[i][-1]))

for paar in werte:
    print(f"({round(paar[0], 5)}, {round(paar[1], 5)})")

integral = 0
for i in range(len(werte)-1):
    integral += (werte[i+1][1]+werte[i][1])/2 * (werte[i+1][0]-werte[i][0])

print(f"Integral = {round(integral, 5)}")