"""Lorenz demo endpoint check. Python 3, standard library only."""
from math import dist, isfinite

def f(p):
    x, y, z = p
    return (10*(y-x), x*(28-z)-y, x*y-(8/3)*z)

def add(p, q, h):
    return tuple(x+h*y for x, y in zip(p, q))

def solve(p, h):
    for _ in range(round(30/h)):
        a = f(p)
        b = f(add(p, a, h/2))
        c = f(add(p, b, h/2))
        d = f(add(p, c, h))
        p = tuple(p[i]+h*(a[i]+2*b[i]+2*c[i]+d[i])/6 for i in range(3))
    assert all(map(isfinite, p))
    return p

for h in (0.00025, 0.000125):
    a = solve((1, 1, 1), h)
    b = solve((1.000001, 1, 1), h)
    print(f"dt={h}: separation at t=30 = {dist(a, b):.6f}")