function dormandPrince(f, t0, y0, t1, h, rtol, atol, maxIter = 10000) {
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 = [35 / 384, 0, 500 / 1113, 125 / 192, -2187 / 6784, 11 / 84, 0];
const dc = [
5179 / 57600,
0,
7571 / 16695,
393 / 640,
-92097 / 339200,
187 / 2100,
1 / 40,
];
let t = t0;
let y = y0;
let iter = 0;
const ys = [y0];
let k = []
k.push(f(t,y))
for(i=1; i<7; i++){
k.push(math.multiply(a[i-1], math.multiply(k, f(t+a[i-1]))))
}
while (t < t1 && iter < maxIter) {
const k = Array(7).fill(0).map((_, i) => {
if (i === 0) return f(t, y);
const sum = a[i - 1].reduce((acc, coef, j) => math.add(acc, math.multiply(coef, k[j])), 0);
return f(t + a[i - 1][0] * h, math.add(y, math.multiply(h, sum)));
});
const yNext = math.add(y, math.multiply(h, math.dot(c, k)));
const yErr = math.multiply(h, math.dot(dc, k));
const errNorm = math.norm(yErr) / (math.norm(yNext) * rtol + atol);
if (errNorm <= 1) {
t += h;
y = yNext;
ys.push(y);
}
h *= 0.9 * Math.pow(errNorm, -1 / 5);
iter += 1;
}
if (iter === maxIter) {
console.warn('Maximum number of iterations reached.');
}
return ys;
}