The dsolve command with the options numeric and method=mebdfi finds a numerical solution of a DAE system using the Modified Extended Backward Differentiation Equation Implicit method. This is a stiff method, so can handle stiff problems efficiently, but should only be used for DAE of index 2 or lower. It can be used to solve regular ODE problems, but use of regular ODE solvers is recommended for that purpose (see dsolve[numeric,IVP]).

The following options are available for the mebdfi method:

'output' and 'known'

The 'output' option specifies the output from dsolve, and the known option specifies user-defined known functions. These options are discussed in dsolve[numeric].

'startinit' and 'optimize'

These options control the method and behavior of the computation, the 'startinit' option is the same as for IVP problems, and is discussed in dsolve[numeric,IVP], while the 'optimize' option is discussed in dsolve[numeric].

'maxfun'

This option specifies the maximum number of steps taken to obtain the requested solution. For this direct solver, this is the closest parallel to dsolve[maxfun]. By default this option is disabled.

'abserr' and 'relerr'

These options specify the desired accuracy of the solution, and are discussed in dsolve[Error_Control]. The default values for mebdfi are $\mathrm{abserr}=1.\times {10}^{\mathrm{-7}}$ and $\mathrm{relerr}=1.\times {10}^{\mathrm{-6}}$.

'minstep', 'maxstep', and 'initstep'

These options provide finer control over the step size used in the method, and are also discussed in dsolve[Error_Control]. By default $\mathrm{minstep}$ is determined within the computation, $\mathrm{maxstep}$ is disabled, and $\mathrm{initstep}=\mathrm{relerr}$. When setting $\mathrm{minstep}$ ensure that the value is adequately small to allow for the first few steps, which are computed at a lower order.

'maxord'

This option specifies the maximum order of the method used in the course of the computation as an integer between 2 and 8. The default value is $8$, but for difficult problems it may be necessary to specify this as $4$ or $5$. Lower order gives better stability, but is less efficient.

The computation may return with an error message corresponding to an error condition of the mebdfi procedure.

Results can be plotted using the function odeplot in the plots package.

$\mathrm{dsys}≔\left\{\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right)\right)\mathrm{diff}\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right)\,t\right)+\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right)\right)\mathrm{diff}\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right)\,t\right)\,\mathrm{x1}\left(t\right)\mathrm{diff}\left(\mathrm{x1}\left(t\right)\,t\right)+\mathrm{y1}\left(t\right)\mathrm{diff}\left(\mathrm{y1}\left(t\right)\,t\right)\,\mathrm{diff}\left(\mathrm{y1}\left(t\right)\,t\,t\right)+9.8+2\mathrm{\λ1}\left(t\right)\mathrm{y1}\left(t\right)+2\mathrm{\λ2}\left(t\right)\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right)\right)\,\mathrm{diff}\left(\mathrm{x1}\left(t\right)\,t\,t\right)+2\mathrm{\λ1}\left(t\right)\mathrm{x1}\left(t\right)+2\mathrm{\λ2}\left(t\right)\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right)\right)\,\mathrm{diff}\left(\mathrm{y2}\left(t\right)\,t\,t\right)+9.8-2\mathrm{\λ2}\left(t\right)\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right)\right)\,\mathrm{diff}\left(\mathrm{x2}\left(t\right)\,t\,t\right)-2\mathrm{\λ2}\left(t\right)\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right)\right)\right\}$

$\mathrm{dsys}≔\left\{\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right)\right)\left(\frac{\ⅆ}{\ⅆt}\mathrm{x1}\left(t\right)-\frac{\ⅆ}{\ⅆt}\mathrm{x2}\left(t\right)\right)+\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right)\right)\left(\frac{\ⅆ}{\ⅆt}\mathrm{y1}\left(t\right)-\frac{\ⅆ}{\ⅆt}\mathrm{y2}\left(t\right)\right)\,\mathrm{x1}\left(t\right)\left(\frac{\ⅆ}{\ⅆt}\mathrm{x1}\left(t\right)\right)+\mathrm{y1}\left(t\right)\left(\frac{\ⅆ}{\ⅆt}\mathrm{y1}\left(t\right)\right)\,\frac{{\ⅆ}^2}{\ⅆt^2}\mathrm{x2}\left(t\right)-2\mathrm{\λ2}\left(t\right)\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right)\right)\,\frac{{\ⅆ}^2}{\ⅆt^2}\mathrm{x1}\left(t\right)+2\mathrm{\λ1}\left(t\right)\mathrm{x1}\left(t\right)+2\mathrm{\λ2}\left(t\right)\left(\mathrm{x1}\left(t\right)-\mathrm{x2}\left(t\right)\right)\,\frac{{\ⅆ}^2}{\ⅆt^2}\mathrm{y2}\left(t\right)+9.8-2\mathrm{\λ2}\left(t\right)\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right)\right)\,\frac{{\ⅆ}^2}{\ⅆt^2}\mathrm{y1}\left(t\right)+9.8+2\mathrm{\λ1}\left(t\right)\mathrm{y1}\left(t\right)+2\mathrm{\λ2}\left(t\right)\left(\mathrm{y1}\left(t\right)-\mathrm{y2}\left(t\right)\right)\right\}$

$\mathrm{ics}≔\left\{\mathrm{x1}\left(0\right)=0\,\mathrm{x2}\left(0\right)=0\,\mathrm{y1}\left(0\right)=-1\,\mathrm{y2}\left(0\right)=-\frac{3}{2}\,\mathrm{D}\left(\mathrm{x1}\right)\left(0\right)=-3\,\mathrm{D}\left(\mathrm{x2}\right)\left(0\right)=4\,\mathrm{D}\left(\mathrm{y1}\right)\left(0\right)=0\,\mathrm{D}\left(\mathrm{y2}\right)\left(0\right)=0\right\}$

$\mathrm{ics}≔\left\{\mathrm{x1}\left(0\right)=0\,\mathrm{x2}\left(0\right)=0\,\mathrm{y1}\left(0\right)=\mathrm{-1}\,\mathrm{y2}\left(0\right)=-\frac{3}{2}\,\mathrm{D}\left(\mathrm{x1}\right)\left(0\right)=\mathrm{-3}\,\mathrm{D}\left(\mathrm{x2}\right)\left(0\right)=4\,\mathrm{D}\left(\mathrm{y1}\right)\left(0\right)=0\,\mathrm{D}\left(\mathrm{y2}\right)\left(0\right)=0\right\}$

$\mathrm{dsol1}≔\mathrm{dsolve}\left(\mathrm{dsys}\mathbf{union}\mathrm{ics}\,\mathrm{numeric}\,\mathrm{method}=\mathrm{mebdfi}\right)$

$\mathrm{dsol1}≔\mathbf{proc}\left(\mathrm{x\_mebdfi}\right)...\mathbf{end\; proc}$

$\mathrm{t1}≔\mathrm{time}\left(\right)\:$

$\mathrm{dsol1}\left(10\right)$

$\left[t=10.\,\mathrm{\λ1}\left(t\right)=8.05374866836774572\,\mathrm{\λ2}\left(t\right)=\mathrm{-6.40533775902082869}\,\mathrm{x1}\left(t\right)=\mathrm{-0.355531961724009915}\,\frac{\ⅆ}{\ⅆt}\mathrm{x1}\left(t\right)=\mathrm{-1.25978921249117670}\,\mathrm{x2}\left(t\right)=\mathrm{-0.492864788890956418}\,\frac{\ⅆ}{\ⅆt}\mathrm{x2}\left(t\right)=\mathrm{-1.05303335398026632}\,\mathrm{y1}\left(t\right)=\mathrm{-0.934664111102251338}\,\frac{\ⅆ}{\ⅆt}\mathrm{y1}\left(t\right)=0.479204587580475905\,\mathrm{y2}\left(t\right)=\mathrm{-0.453894585514793392}\,\frac{\ⅆ}{\ⅆt}\mathrm{y2}\left(t\right)=0.538264833117563346\right]$

$\mathrm{time}\left(\right)-\mathrm{t1}$

$0.749$

$\mathrm{dsol2}≔\mathrm{dsolve}\left(\mathrm{dsys}\mathbf{union}\mathrm{ics}\,\mathrm{numeric}\,\mathrm{method}=\mathrm{mebdfi}\,\mathrm{maxord}=4\right)$

$\mathrm{dsol2}≔\mathbf{proc}\left(\mathrm{x\_mebdfi}\right)...\mathbf{end\; proc}$

$\mathrm{t2}≔\mathrm{time}\left(\right)\:$

$\mathrm{dsol2}\left(10\right)$

$\left[t=10.\,\mathrm{\λ1}\left(t\right)=8.05124380921260041\,\mathrm{\λ2}\left(t\right)=\mathrm{-6.42065568374894546}\,\mathrm{x1}\left(t\right)=\mathrm{-0.357252010506513173}\,\frac{\ⅆ}{\ⅆt}\mathrm{x1}\left(t\right)=\mathrm{-1.25072569463780869}\,\mathrm{x2}\left(t\right)=\mathrm{-0.492112228851291800}\,\frac{\ⅆ}{\ⅆt}\mathrm{x2}\left(t\right)=\mathrm{-1.04698380844525074}\,\mathrm{y1}\left(t\right)=\mathrm{-0.934008348789023168}\,\frac{\ⅆ}{\ⅆt}\mathrm{y1}\left(t\right)=0.478394293227407508\,\mathrm{y2}\left(t\right)=\mathrm{-0.452537634291989843}\,\frac{\ⅆ}{\ⅆt}\mathrm{y2}\left(t\right)=0.535462511444270373\right]$

$\mathrm{time}\left(\right)-\mathrm{t2}$

$1.234$

Cash, J.R. "The Integration of stiff IVP in ODE using modified extended BDF." Computers and Mathematics with Applications. Vol. 9. (1983): 645-657.

Cash, J.R., and Considine, S. "An MEBDF code for stiff IVP." ACM Trans Math Software. 1992: 142-158.