def blackbox(f, t0=0.0, y0=1.0, nsteps=100, h=1e-2):
"""
Method to solve an IVP of shape
y' = f(t,y); y(t_0) = y_0
where len(y_0) determines the problem dimension.
Arguments:
f (callable) Right hand side of ODE in standard form.
t_0 (float) Initial time value. Default is 0.
y_0 (np.array) Initial value. Default is 1.
nsteps (int) Number of integration steps. Default is 100.
h (float) Step length. Default is 1e-2.
Returns:
timepoints ([float]) Computed time steps.
yapprox ([np.array]) Computed solution approximations.
"""
# init output variables
timepoints = [t0]
yapprox = [y0]
# init iteration variables
t = t0
y = y0
for i in range(nsteps):
# intermediate computations for current step
k1 = f(t, y)
k2 = f(t + (1/5)*h, y + h*((1/5)*k1))
k3 = f(t + (3/10)*h, y + h*((3/40)*k1 + (9/40)*k2))
k4 = f(t + (4/5)*h, y + h*((44/45)*k1 - (56/15)*k2 + (32/9)*k3))
k5 = f(t + (8/9)*h, y + h*((19372/6561)*k1 - (25360/2187)*k2 + (64448/6561)*k3 - (212/729)*k4))
k6 = f(t + h, y + h*((9017/3168)*k1 - (355/33)*k2 + (46732/5247)*k3 + (49/176)*k4 - (5103/18656)*k5))
# perform step
y = y + h*((35/384)*k1 + (500/1113)*k3 + (125/192)*k4 - (2187/6784)*k5 + (11/84)*k6)
t = t + h
timepoints.append(t)
yapprox.append(y)
return timepoints, yapprox