gistfile1.txt

 let t0 = performance.now();
    let N = 3000; //iteracje
    let TAU = 0.01; //TAU

    let particles = 10;

    let damping = 0.5;

    function RK4(tau, n) {
        let t0 = 0.0;
        let y0 = [];

        for (let z = 0; z < 4 * particles; z++) {
            let random = Math.random() * 8 - 4;
            y0 = [...y0, random]
        }

        //y0 = [1, 0, 0, 0, -1, 0, 0, 0, 0, 0.8, 0, 0];
        // x,y,vx,vy
        let yy = [y0], tt = [t0], i;


        function g(t, y) {
            let g = [];

            for (let c = 0; c < particles; c++) {
                let temp = c * 4;

                let y1, y2, y3, y4;
                y1 = y[2 + temp];
                y2 = y[3 + temp];
                y3 = -y[temp];
                y4 = -y[1 + temp];
                for (let a = 0; a < particles; a++) {
                    if (a === c) {
                        continue;
                    }
                    y3 += ((y[temp] - y[a * 4]) / (Math.pow(Math.pow(y[temp] - y[a * 4], 2) + Math.pow(y[1 + temp] - y[1 + a * 4], 2), 3 / 2)));
                    y4 += ((y[1 + temp] - y[1 + a * 4]) / (Math.pow(Math.pow(y[temp] - y[a * 4], 2) + Math.pow(y[1 + temp] - y[1 + a * 4], 2), 3 / 2)));

                }
                y3 -= damping * y[2 + temp];
                y4 -= damping * y[3 + temp];

                g = [...g, y1, y2, y3, y4];

            }
            return g;
        }


        function K1(t, y) {
            return g(t, y);
        }

        function K2(t, y) {
            let temp = K1(t, y);
            let yy = [];
            for (let c = 0; c < particles; c++) {
                yy[c * 4] = y[(c * 4)] + (tau / 2) * temp[(c * 4)];
                yy[1 + c * 4] = y[1 + c * 4] + (tau / 2) * temp[1 + c * 4];
                yy[2 + c * 4] = y[2 + c * 4] + (tau / 2) * temp[2 + c * 4];
                yy[3 + c * 4] = y[3 + c * 4] + (tau / 2) * temp[3 + c * 4];
            }
            return g(t + tau / 2, yy);
        }

        function K3(t, y) {
            let temp = K2(t, y);
            let yy = [];
            for (let c = 0; c < particles; c++) {
                yy[(c * 4)] = y[(c * 4)] + (tau / 2) * temp[(c * 4)];
                yy[1 + c * 4] = y[1 + c * 4] + (tau / 2) * temp[1 + c * 4];
                yy[2 + c * 4] = y[2 + c * 4] + (tau / 2) * temp[2 + c * 4];
                yy[3 + c * 4] = y[3 + c * 4] + (tau / 2) * temp[3 + c * 4];
            }
            return g(t + tau / 2, yy);
        }

        function K4(t, y) {
            let temp = K3(t, y);
            let yy = [];
            for (let c = 0; c < particles; c++) {
                yy[(c * 4)] = y[(c * 4)] + tau * temp[(c * 4)];
                yy[1 + c * 4] = y[1 + c * 4] + tau * temp[1 + c * 4];
                yy[2 + c * 4] = y[2 + c * 4] + tau * temp[2 + c * 4];
                yy[3 + c * 4] = y[3 + c * 4] + tau * temp[3 + c * 4];
            }
            return g(t + tau, yy);
        }

        for (i = 0; i < n; i++) {
            let yn = [];
            for (let c = 0; c < particles; c++) {
                yn[(c * 4)] = y0[(c * 4)] + (tau / 6) * (K1(t0, y0)[(c * 4)] + 2 * K2(t0, y0)[(c * 4)] + 2 * K3(t0, y0)[(c * 4)] + K4(t0, y0)[(c * 4)]);
                yn[1 + c * 4] = y0[1 + c * 4] + (tau / 6) * (K1(t0, y0)[1 + c * 4] + 2 * K2(t0, y0)[1 + c * 4] + 2 * K3(t0, y0)[1 + c * 4] + K4(t0, y0)[1 + c * 4]);
                yn[2 + c * 4] = y0[2 + c * 4] + (tau / 6) * (K1(t0, y0)[2 + c * 4] + 2 * K2(t0, y0)[2 + c * 4] + 2 * K3(t0, y0)[2 + c * 4] + K4(t0, y0)[2 + c * 4]);
                yn[3 + c * 4] = y0[3 + c * 4] + (tau / 6) * (K1(t0, y0)[3 + c * 4] + 2 * K2(t0, y0)[3 + c * 4] + 2 * K3(t0, y0)[3 + c * 4] + K4(t0, y0)[3 + c * 4]);
            }

            t0 += tau;
            yy.push(yn);
            y0 = yn;
            tt.push(t0);
        }

        return [tt, yy];
    }