---

y_mod23 = np.zeros(n_intervals+1)
y_mod23[0] = y0_23
for i in range(n_intervals):
    f1 = f23(x_grid23[i], y_mod23[i])
    y_mid = y_mod23[i] + (h23/2) * f1
    f_mid = f23(x_grid23[i] + h23/2, y_mid)
    y_mod23[i+1] = y_mod23[i] + h23 * f_mid