My documentation shows that at v2.5.2 I fixed a bug in CombineTrajectories which was leading to extremely slow calculations, and then fixed another bug in the same routine at v2.5.4, and again for v2.6b2. As I recall, I was seeing mysterious overshoots in the trajectory. I traced these to the FunctionInterpolation function of Mathematica, which I found to be unreliable. So I recoded to avoid FunctionInterpolation, but the efficiency of the code seems to have been compromised at the same time. I tend not to use so many time intervals, so I must have missed its poor efficiency in that case.

Unfortunately I am reluctant to open that can of worms again, given my current commitments.

If anyone volunteers to look at this, I have on file the test notebooks that I used to fix the bugs mentioned above.

If I recall, the line mentioned above:

`\[Rho]intervalfuncs = Function[t,Evaluate[Assuming[(First[#] < t < Last[#]), Refine[\[Rho]func[t]]]]]& /@tintervals;`

Has the following purpose.

tintervals contains the set of time intervals for the Piecewise function, i.e. {..,{tb,ta}}

\[Rho]func is the Piecewise function for the density operator Rho, which may be evaluated at any time over the whole trajectory

\[Rho]intervalfuncs is the set of functions for Rho, which apply for the individual intervals, i.e. {..Rhofunc2, Rhofunc1}.

I obtain these from Rhofunc by using Refine and Assuming. This has the benefit of reliability and generality. I found this method to be convenient but there might well be a more efficient method.