Maple and Mathematica disagree using dsolve for system of ODE initial value problem

Question

I'm trying to solve a system of ODEs, symbolically, using Maple and Mathematica. I'm actually comparing the analytic solution with a numerical one (this is for a large simulation project, of which this is a part), and what I find is that the Mathematica solution appears to follow the numerical solutions I'm running, whereas the Maple solution does not. The solutions in this case are too long to be computationally feasible when compared with a numerical solution anyway, so I won't be using either answer, but I'm very curious as to why there is a disagreement. The system of equations is not complex:

$$y_1=A(e^{-at}-e^{-bt})$$ $$\dot{y_2}=a_1y_1-a_2y_2$$ $$\dot{y_3}=a_2y_2-a_3y_3$$

The solution to $y_2$ is the same between Maple and Mathematica, but $y_3$ is very different. I'm aware that sometimes approximate solutions are used, but no warnings came up for this. Is there some other uniqueness in how the two packages handle symbolic solutions that I'm missing?

For completeness, my Mathematica code is: (https://i.stack.imgur.com/3qcr9.png)

And for Maple: (https://i.stack.imgur.com/pjAiS.png)

Any thoughts? Many thanks!

EDIT: While the initial values in general are variables, for purposes of comparison, I set them all to 0.

EDIT 2: Even though the form of the equations are very different, I know they could be mathematically identical. I think they are different, though, because I substituted values for all the parameters and get two different numbers as output.

Answer 1

The problem was one of those RATS (Read All The Screen) errors. Thanks to the comment by Mariusz I was able to find the bug. So this isn't really a good question after all. But I would like to keep this silly question up because of the answer that did come. It seems good to know.

Answer 2

This is not an answer, just an extended comment. I have not look at your codes, but if you want to find an analytic form of the solution it is quite easy to do by hand. The system you are considering is $$ \begin{bmatrix}\dot y_2\\\dot y_3\end{bmatrix}= \begin{bmatrix}-a_2 & 0\\a_2 & -a_3\end{bmatrix}\begin{bmatrix} y_2\\y_3\end{bmatrix}+\begin{bmatrix}a_1\\0\end{bmatrix}y_1. $$ It takes the form of linear equation $$ \dot x=Ax+Bu $$ with the general solution being $$ x(t)=e^{At}x(0)+e^{At}\int_0^te^{-As}Bu(s)\,ds. $$ In your case, $u(s)=y_1(s)$ is a scalar function of exponentials, so the integral can be easily calculated. For example, $$ \int_0^te^{-As}Be^{-as}\,ds=\int_0^te^{-(A+aI)s}\,ds\cdot B=(A+aI)^{-1}(I-e^{-(A+aI)t})\cdot B. $$