import matplotlib.pyplot as plt
import math


def f(x):
    return 1 / math.sqrt(2) * (x - 1) / math.sqrt(x - math.log(x))


def df(x, y):
    return y / (x - 1) - y ** 3 / (x * x - x)


a = 2
b = 20

n = 5

x0 = a
y0 = f(x0)

x1 = []
y1 = []
y2 = []
h1 = (b - a) / n

x1.append(x0)
y1.append(y0)
y2.append(y0)

for i in range(1, n + 1, 1):
    k1 = h1 * df(x1[-1], y1[-1])
    k2 = h1 * df(x1[-1] + h1 / 2, y1[-1] + k1 / 2)
    k3 = h1 * df(x1[-1] + h1 / 2, y1[-1] + k2 / 2)
    k4 = h1 * df(x1[-1] + h1, y1[-1] + k3)

    y1.append(y1[-1] + 1 / 6 * (k1 + 2 * k2 + 2 * k3 + k4))

    x1.append(x0 + i * h1)
    y2.append(f(x1[-1]))


eps_lim = 0.01
eps_max = 0

for i in range(1, n + 1):
  l = x1[i - 1]
  r = x1[i]
  yl = y1[i - 1]
  yr = y1[i]
  ycur = yl
  xcur = l
  while xcur < r:
    xcur += (r - l) / 100
    ycur += (yr - yl) / 100
    eps_max = max(eps_max, abs(f(xcur) - ycur))



print(f"eps = {eps_max}")

print(f"n = {n},\nh = {h1}\ny(b) = {y1[-1]}")
plt.plot(x1, y1, color='red')
plt.plot(x1, y2, color='green')
plt.show()