function ode45(f, T, y0){
const tol = 1e-4
const hmax = 1
const t0 = T[0]
const tf = T[1]
const max_iter = 10_000
const minSteps = 8
let h = T[2] ? T[2] : math.divide(math.subtract(T[1],T[0]), minSteps)
let t = [t0]
let y = [y0]
const a = [
[],
[1/5],
[3/40, 9/40],
[44/45, -56/15, 32/9],
[19372/6561, -25360/2187, 64448/6561, -212/729],
[9017/3168, -355/33, 46732/5247, 49/176, -5103/18656],
[35/384, 0, 500/1113, 125/192, -2187/6784, 11/84]
]
const c = [0, 1/5, 3/10, 4/5, 8/9, 1, 1]
const b = [35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0]
const bp = [5179/57600, 0, 7571/16695, 393/640, -92097/339200, 187/2100, 1/40]
const delta_b = math.subtract(b, bp)
let n = 0
let iter = 0
while (notFinished(t[n], tf, h) && (iter < max_iter)){
let k = []
h = trimStep(t[n], tf, h)
k.push(f(t[n], y[n]))
k.push(f(math.add(t[n], math.dotMultiply(c[1], h)),
math.add(
y[n],
math.dotMultiply(
h,
math.dotMultiply(a[1][0], k[0])
)
)
))
k.push(f(
math.add(t[n], math.dotMultiply(c[2], h)),
math.add(
y[n],
math.dotMultiply(
h,
math.add(
math.dotMultiply(a[2][0], k[0]),
math.dotMultiply(a[2][1], k[1])
)
)
)
))
k.push(f(
math.add(t[n], math.dotMultiply(c[3], h))
,
math.add(
y[n],
math.dotMultiply(
h,
math.add(
math.dotMultiply(a[3][0], k[0]),
math.dotMultiply(a[3][1], k[1]),
math.dotMultiply(a[3][2], k[2])
)
)
)
))
k.push(f(
math.add(t[n], math.dotMultiply(c[4], h)),
math.add(
y[n],
math.dotMultiply(
h,
math.add(
math.dotMultiply(a[4][0], k[0]),
math.dotMultiply(a[4][1], k[1]),
math.dotMultiply(a[4][2], k[2]),
math.dotMultiply(a[4][3], k[3])
)
)
)
))
k.push(f(
math.add(t[n], h),
math.add(
y[n],
math.dotMultiply(
h,
math.add(
math.dotMultiply(a[5][0], k[0]),
math.dotMultiply(a[5][1], k[1]),
math.dotMultiply(a[5][2], k[2]),
math.dotMultiply(a[5][3], k[3]),
math.dotMultiply(a[5][4], k[4])
)
)
)
))
k.push(f(
math.add(t[n], h),
math.add(
y[n],
math.dotMultiply(
h,
math.add(
math.dotMultiply(a[6][0], k[0]),
math.dotMultiply(a[6][1], k[1]),
math.dotMultiply(a[6][2], k[2]),
math.dotMultiply(a[6][3], k[3]),
math.dotMultiply(a[6][4], k[4]),
math.dotMultiply(a[6][5], k[5])
)
)
)
))
const TE = math.abs(
math.add(
math.dotMultiply(delta_b[0], k[0]),
math.dotMultiply(delta_b[2], k[2]),
math.dotMultiply(delta_b[3], k[3]),
math.dotMultiply(delta_b[4], k[4]),
math.dotMultiply(delta_b[5], k[5]),
math.dotMultiply(delta_b[6], k[6])
)
)
const TEnum = math.isUnit(TE) ? TE.toNumber(TE.formatUnits()) : TE
const delta =
0.84 *
math.dotPow(
math.dotDivide(tol, TEnum),
(1/5)
)
if ( TEnum < tol ) {
t.push(math.add(t[n], h))
y.push(
math.add(
y[n],
math.dotMultiply(
h,
math.multiply(b, k)
)
)
)
n += 1
}
if (delta <= 0.1){
h = math.multiply(h, 0.1)
}
else if( delta >= 4)
{
h = math.multiply(h, 4)
}
else{
h = math.multiply(h, delta)
}
if ( h > hmax){
h = hmax
}
//h = math.multiply(h, delta)
iter ++
}
return {t,y}
}