I took some time to make it work (or at least address some major issues) in an online compiler. Here are my findings:

Scale Issues

This was an easy problem to fix. As i already mentioned in my comment, these formulas only work in the International System of Units (SI units), and when you scale down meters for your simulation (to avoid getting huge numbers in your rendering logic i assume, or to make them easier to read), you would also have to scale down everything else.

Because many formulas are not linear (for example: the gravity experienced by an object is depending on the square of the distance [1], so if you half the simulation distance and half your simulation mass, the result doesn't match anymore.

Therefore, i'd strongly recommend against scaling at all, at least for the physics.

In my project, i have seperated rendering and physics. You can define a factor (like 1 AU or ~10^-12, depending on your needs). For quick testing, i defined a 1 AU factor and applied it to your data (and my made up data):

const astronomicalUnit = 1.496 * (10**11);
const bodies = [
    {
        position: [0, 0.8 * astronomicalUnit],
        velocity: [0, 0],
        force: [0, 0],
        mass: 1.989 * 10 ** 30,
        radius: 3.5,
        trailPositions: [],
        colour: [1, 1, 0.8, 1],
        parentIndex: -1,
        name: "Sun",
    },
    {
        position: [29.7 * astronomicalUnit, 0], // Approximate distance at lowest point in orbit
        velocity: [0, 6.1 * 10 ** 3], // Orbital speed in m/s (approximate data for lowest point aswell)
        force: [0, 0],
        mass: 1.309 * 10 ** 22, // kg
        radius: 0.0064, // Scaled radius for visualization; actual radius ~1.1883 × 10^6 m
        trailPositions: [],
        colour: [0.6, 0.7, 0.8, 1], // Pale bluish-grey
        parentIndex: 0, // Assuming Sun is index 0
        name: "Pluto", // I picked pluto because it has well known orbital data and an eccentricity of 0.25, 
                       // which should make the second focus visually distinct
    }
];

And after all calculations are done, you can divide out the same factor to get a more readable and renderable result. It also enables you to use more even more real-world constants:

const gravity = 6.674 * (10**(-11)); // real world gravitational constant
findSecondFocus(1);

// Within findSecondFocus:
console.log("Semi-major axis:", (a / astronomicalUnit)); // Instead of printing a directly, for example

This already fixes the calculation of the semi-major axis!

To summarize: use realistic values, if you want realistic results (alternatively: experiment to find consistent values for an alternate universe, but that will take time and disable you from just looking up data). Most relevant for your project: meters, kilograms and seconds.

Implementation Issues

Here:

    // The eccentricity vector formula is: e = (v × h)/μ - r/|r|
    const rvDot = relativeSpatiumVector[0] * relativeVelocityVector[0] +
              relativeSpatiumVector[1] * relativeVelocityVector[1];

You write cross product in your comment, but use the dot product. You also use the dot product to calculate h, the angular momentum vector. Unfortunately, it takes quite a bit of effort to fix this one.

The cross product of two vectors produces a vector perpendicular to both input vectors [2]. Where does it go for 2D vectors? Outside of your plane of simulation.

Thats quite unfortunate, but we can cheese our way around.

First, i made some helpers for both 2D and 3D cross products:

// Seperate definition of a cross-product helper, so code is easier to read
function cross2D(a, b) {
    return a[0] * b[1] - a[1] * b[0];
}
function cross3D(a, b) {
    return [
        a[1] * b[2] - a[2] * b[1],
        a[2] * b[0] - a[0] * b[2],
        a[0] * b[1] - a[1] * b[0]
    ];
}

Then, i replaced the code for eccentricity vector calculation, i'll explain afterwards:

    // The eccentricity vector formula is: e = (v × h)/μ - r/|r|
    const rUnit = [
        relativeSpatiumVector[0] / r,
        relativeSpatiumVector[1] / r
    ];
    
    const angular_z = cross2D(relativeSpatiumVector, relativeVelocityVector);
    const angularMomentumVector = [0,0,angular_z]; // This is the "h"
    const liftedVelocityVector = [relativeVelocityVector[0], relativeVelocityVector[1], 0];
    const vxh = cross3D(liftedVelocityVector, angularMomentumVector);
    const eccentricityVector = [
        vxh[0] / mu_sim - rUnit[0],
        vxh[1] / mu_sim - rUnit[1],
    ]; // (v × h)/μ - r/|r|

Your rUnit looked fine, so i reused it. I created a angular velocity 3D vector angularMomentumVector by assuming everything on the 2D plane to be zero, which i can do because it has to be perpendicular to two vectors on this plane.

Then, we need to get the velocity into 3D (liftedVelocityVector) aswell. Thats easy, because it just doesn't move in the z direction.

Then, we get the cross product in vxh, and can finally apply the formula you already had in your comment.

We can ignore the z component (vxh[2]), because the cross product must be perpendicular to the angularMomentumVector, which only has z components.

Everything else in your code was perfectly fine, so well done!

Test results

With the data from earlier in the answer and these updated console logs:

    console.log("Second Focus coordinates:", secondFocus[0] / astronomicalUnit, ", ", secondFocus[1] / astronomicalUnit);
    console.log("Eccentricity:", eccentricityScalar);
    console.log("Semi-major axis:", (a / astronomicalUnit));

I get these results:

Second Focus coordinates: -19.369704292780035 ,  -1.321742876573199
Eccentricity: 0.2472841556295451
Semi-major axis: 39.39913738651615

Compared to Wikipedia Data, thats ~0.0015 off in eccentricity, and ~0.083 AU off in the semi-major axis. I blame the inaccuracy on my rounded input data and the fact we clipped off its entire inclination.

I could not find a reference value for the second focus, but it seems plausible.

Thanks for the fun challange and good look with your project!