How to choose initial state vectors for a resonant orbit numerical simulation?
I don't know if there are rules of thumb for initializing resonant orbital simulations.  I might suggest:

Setting initial state vectors' Z-Components to zero.  Keeping the whole simulation in the X-Y plane can help with debugging and visualization.

Choosing initial state vectors that are equivalent to orbits with non-zero eccentricity. Orbital resonances tend to amplify or dampen eccentricities, but the effects can be multiplicative. So if you choose circular orbits, you may see less of an effect.

Choosing initial state vectors that line up apogees and perigees with closest contacts.

Running multiple simulations with slight perturbations to initial state vectors to test stability.

Some interesting initial experiments might be to see if you could replicate the instability of the L1 and L2 points of an orbit.  Or perhaps put some little masses in the Kirkwood Gaps and see if they get expelled.


Notes:

I tried writing some of my own numerical simulations to simulate natural planetary flybys. I used an algorithm outlined by Voesenek, but I was getting round-off errors that shot my simulation. Instead of investing more time to fix my simulation, I ended up using Universe Sandbox 2, that I paid for and downloaded from Steam.

Universe Sandbox probably doesn't have the fidelity you need for resonance stability simulations.  Gallardo wrote a pretty slick long term orbital integration software package called ORBE, that is free to the public for educational purposes.  I wish I had read Tabaré Gallardo's 2017 Exploring the orbital evolution of planetary systems before writing my own simulator, since he uses a neat trick by Encke to account for the Sun's gravity analytically via 2-body problem equations while separately treating the inter-planetary gravitational accelerations numerically. This allows larger time steps and prevents some kinds of round-off error.  Gallardo's software is designed to test long term stellar system stability.  Also, it uses orbital elements rather than state vectors, which are a bit less intuitive to me. Gallardo's used of ORBE to run Kirkwood Gap stability simulations is totally AWESOME!!  Here is a capture of one of the figures from Gallardo's 2017: