The Problem
I have a differential equation of the following form $$\dot{y} = A(y)y+B(y)$$ Provided that $A$ does not depend on $y$ (or depends very weakly on $y$), I can approximate by using the following form: $$y(t)=e^{-At}y_0+\int^t_0e^{-(t-\tau)A}B(y)d\tau$$ Re-introducing the dependence of $A$ on $y$, the above can be discretised as follows:$$y_{k+1}=e^{-A(y_k)\Delta t}y_k + \Delta t\ \varphi_1(-\Delta t A)B(y_k)$$ where $\varphi_1(z) = \frac{e^z-z}{z}$ and $\Delta t$ is the time step.
I can confirm that the above approximation is accurate to make, since the result has been compared with another typical implicit solver, and the difference between the two is acceptable (nearly negligible).
The Question
However, my question is, can we get an expression for the local truncation error of the above? I have looked in this paper, which extensively describes such type of integrators, however I can't seem to be able to understand how to obtain such an expression.
My thinking is that if I can estimate the local truncation error of the above, I can find a way to automatically adjust the timestep in accordance to my desired tolerances.