The only documentation I know about for the implementation of
ode23t is in the paper which documents the implementation of ode23tb, the TRBDF2 method in MATLAB.
As usually implemented, the trapezoidal rule is not strictly a one-step
method because the truncation error estimate makes use of two previously computed solution values. The method is not efficient for very stiff problems because it is not strongly stable.
and it cites The SPICE book. Two things here. First of all, the Trapezoidal rule, also known as Crank-Nicholson, is not L-stable and is generally not suited for highly stiff equations (but this lack of dampening can give undesirable outcomes in some models). This is exacerbated by the instability of the error estimate, which is from a truncation error estimate. Basically, you save some previous points and use derivative approximations in order to match the leading truncation error terms, and then use this as the input to standard P or PI adaptive algorithms.
For an open source (MIT licensed) implementation of this algorithm, you can take a look at the
Trapezoid method from DifferentialEquations.jl. The implementation can be found in the source code here and it was derived from the description in the SPICE book reference given by Shampine. Note that due to licensing issues this has not been directly checked against the MATLAB source code, only directly against these sources given by Shampine.
Responding to Edits
Continuing the conversation is kind of frowned upon in SO, but I'll do it since it makes sense. For part one, the linked method isn't implicit Euler in the first step. Instead, the adaptivity doesn't exist in the first step since it uses two steps to get a second derivative approximation. However, being a 1-step method, the Trapezoidal rule is still used even when it takes the non-adaptive step. This keeps it always second order (though you can have a single first order step and still be second order).
But... as stated above, MatLab estimates the error by differentiating a cubic interpolant.
Where is it stated that MATLAB does this? MATLAB uses a cubic interpolant error for
dde23, but I have never seen it stated that it does this on
ode23t. In the comments I show the paper trail starting from MATLAB's documentation to Shampine's paper and to the SPICE book. In the SPICE book they derive the truncation error estimate which uses two steps to estimate the third derivative term which is the leading truncation error term. To the best of my knowledge, and not looking at the MATLAB source since I am an open source developer and do not want licensing issues, that is how it's done according to the sources given. If someone has the time, they can double change and confirm or deny this by looking at the MATLAB source.
To answer the question of how you could come up with an error estimate using the interpolant, let me show you that Hermite won't work here. I am just going to point to the live code since I know it's correct and it's easier than TeXing it up. The Hermite interpolant is defined here:
with ($y_0$,k) and ($y_1$,k) as the values at t and t+dt respectively, calculating at the point t+$\theta$dt. To get the derivative of the interpolant, you differentiate by $\theta$ and get
Notice that when $\theta=1$, the derivative approximated by the cubic Hermite polynomial is $k$, which is the derivative estimate at $t+dt$ calculated from $f$ itself. So it doesn't have a natural error estimate at the end points of the interval against the derivative since it's made to perfectly hit those derivative values. When this estimate is used it's called a residual correction method, discussed here, where this derivative estimate is used against extra $f$ evaluations in the middle of the interval to give a continuous derivative error estimate. By using specific points you can estimate the maximum difference, and this done on RK4 is how
ddesd works. RK4 is simple, and an open source implementation of the two-point maximal residual error estimate for it can be found here. But while MATLAB does use this in a different integrator, Shampine writes about it, and you technically can apply this to the Trapezoid method, from what I can tell this is not what's done in