If I haven't made some error, it looks like the geodesic equations for the dss metric should boil down to
$$\dot{t} = \frac{E}{f(r)} \quad \dot{r} = \sqrt{E^2-f(r)}$$
It'd take more work to write ##\dot{t}## and ##\dot{r}## as a function of ##\tau##, however, so it's difficult to directly confirm that these are the correct solutions as-is. I imagine one could use the chain rule, but I haven't done this.
Here E is some constant, representing the energy, and f(r) is, as previously
$$f(r) = 1 - \frac{2a}{r} - br^2$$
It looks like the motion is best understood by the effective potential technique, as used by MTW and on the forurmilab website for the Schwarzschild case, and that the effective potential ##V^2(r)## is just f(r). When b=0, this matches the Schwarzschild effective potential ##1-2a/r## from the fourmilab site / MTW as it should.
The forumilab website is
https://www.fourmilab.ch/gravitation/orbits/
a and b are somewhat inconvenient parameter, I wound up normalizing the event horizon to occur at at r=1, and the cosmologcal horizon at r=C, where C is some constant, and then solves for a and b.
Plotting the effective potential for C=10, I get something that looks like this.
If the object can reach the peak of the effective potential at ##r=\sqrt[3]{55} \approx 3.8##, with a coordinate velocity ##dr/d\tau## greater than 0, it will reach the cosmological horizon (and continue on to infinity). I've stop the graph at the cosmological horizon as the coordinates are too confusing beyond it.