Well it’s really hard to help if there’s no code to compare. But the #1 issue always tends to be that the translation is not 1-1. If you run all of the functions with random (u,p,t) do you get the same exact answer from MATLAB and Julia? In almost all cases like this, I’ve found that to be false, so that’s the first thing I would look into detail with.
Other than that ode23t is very different from solve of the solvers you’re discussing here. Did you try the Trapezoidal counterpart? L-stability is good for “most” equations, but not ones for which the dampening is not beneficial. In which case, the reversibility (and “close to symplecticness”) can be beneficial. So it’s good to try methods that are similar to what you’ve seen work well before.
Lastly, you may want to watch this video which highlights some of the most common issues people have when translating, particularly from MATLAB:
