The key is to use a Repeat Zone (hereafter "RZ"). The math and geometry changes to be repeated N times with time steps of Δt/N must be put inside the RZ.
This is my first time using an RZ, so the details and mis-steps are fresh in my mind. I may have some finer points incorrect. This is StackExchange, where anyone can make corrections or clarifications.
Here is what I did, illustrated.
I had been experimenting with Simulation Zone, following a [tutorial][sztut] on making particles follow the basic physics of motion. At 24fps, the simulation was okay by quickie experimental art standards. But for high quality results in scenes with more complex force fields and interactions, sub-steps would be necessary.
My Geometry Nodes setup was like this:

The Simulation takes in positions, velocities and force fields (just gravity for now), calculates updates to them for the next frame, and updates the geometry. Position is updated using Set Position and velocity, and velocity is updated using gravity (acceleration). Doing this in the most naive way possible, frame by frame, is called Euler integration. The Euler node is a Node Group I wrote earlier.
Using shift-A / Utilities I added a Repeat, appearing as the brown area.
How does this work? Starting values are given to the input sockets on the left side, including geometry and any number of other values you may drag in. The list matches the final outputs on the right side, except there's also an Iterations input, for how many times to repeat the computation inside. For brevity I'll use "N" for the number of iterations.
The starting values are given to the nodes inside (none at the moment) which compute new values. These would be sent to the outputs on the right side, except just not yet -- they are sent back to the inputs for recalculation. Only after N iterations do the results appear to the outside world.
In my case, I want to repeat the Euler Integration node N times, and for the Δposition to be added to the geometry's position for each particle. In concept, I want to put the Euler node and the Set Position node into the RZ. But you can't just drag these into the RZ and let Blender figure out the new connections. (Maybe in a future version?)
It's vital to think about exactly which nodes must be repeated for sub-stepping, and which should be executed only once per frame. In my example, every frame creates new particles to be animated, and removes particles that have aged out of existence. It would be a disaster if these operations were repeated with each sub-step. These operations stay outside of the RZ.
You must manually rework the connections. Begin by drawing a connection from the output of the node feeding Set Position (happens to be Delete Geometry in this example) to the Geometry input of the RZ.

Then connect the RZ's Geometry output to the input of the node following Set Position, which happens to be the Simulation's output.
Set Position now has no geometry connections. Give it new connections to the inner sockets of the Repeat nodes. The flow of Geometry data now goes: Delete node, Repeat left-side node, Set Position, Repeat right-side node, then Simulation output.
Do likewise for the Euler node. During the operation the brown and purple polygons marking the Simulation and Repeat Zones may disappear. They'll re-appear when all the connections are complete. Blender's Geometry Nodes Editor will automatically redraw these as it sees suitable, and will automatically put the Euler and Set Position nodes inside the RZ.
The fully wired-up node setup, with some minor mistakes described later, looks like

Note that it's okay to have data flowing into the RZ without passing through the left-hand inputs node. These are constants, at least for going from one frame to the next, not varying as sub-steps are computed. Time should vary, and if I weren't using a simple constant uniform gravitational field, it would need to be recomputed which would require a force field computations inside the RZ.
The one "time" input in my Euler node is misnamed and should be Delta Time or perhaps Time Step. (Oops!) Since the sub-step needs to use a fraction Δt/N of the Δt between frames, I add a Math node outside the RZ to computer this. The N used here must match the Iterations input of the RZ. To make this easy to change, I created an integer constant node to set N (default 1, never allow to be less than 1).

To my disappointment, things didn't work. I was naming sockets however I liked as I went along. I had two vector sockets at both ends of the RZ which didn't make sense. Just need one for velocity. Once I figured out that they must match between the inputs and outputs of the RZ, and the inputs and outputs of the Simulation Zone, this was easy to fix.

References:
[sztut]: About the simulation, no sub-steps: "Simulate Physics in Blender 3.6 / Particles Acceleration and Collisions / Geometry Nodes Tutorial" by Blender Path on YouTube, https://www.youtube.com/watch?v=7Q6KjOCOALk
For use of RZ to create static geometry (not simulation sub-steps) see "Repeat Zone in Geometry nodes - Blender tutorial ENG" by Andrea Ciani, https://www.youtube.com/watch?v=bF4Y6T0r3hc
Another good tutorial on using RZ for static geometry: What is the Repeat Zone in Blender Geometry Nodes? by Stephen Woods, https://www.youtube.com/watch?v=4C9msC_VvRQ