import numpy as np
from scipy.integrate import solve_ivp


decay_rate = 0.4
initial_value = 2.0
sample_times = np.array([0.0, 1.0, 2.5, 5.0, 10.0])


def exponential_decay(t, y):
    return [-decay_rate * y[0]]


solution = solve_ivp(
    exponential_decay,
    t_span=(sample_times[0], sample_times[-1]),
    y0=[initial_value],
    t_eval=sample_times,
    method="RK45",
    rtol=1e-7,
    atol=1e-10,
)

expected = initial_value * np.exp(-decay_rate * sample_times)
max_abs_error = np.max(np.abs(solution.y[0] - expected))

print(f"success: {solution.success}")
print(f"message: {solution.message}")
print(f"t: {np.array2string(solution.t, precision=2)}")
print(f"y: {np.array2string(solution.y[0], precision=8)}")
print(f"expected: {np.array2string(expected, precision=8)}")
print(f"max_abs_error: {max_abs_error:.2e}")
print(f"function_calls: {solution.nfev}")