# Does Mathematica have a command analogous to ode45 of MATLAB?

Does anybody know if Mathematica has an analogue of MATLAB's ode45 command? I need to solve a second order coupled ODE system of equations.

• why do you need a specific integrator? The automatic method selection of Mathematica does this for you in a very sophisticated manner for a much wider range of DEs where it is not easy to know in advance which combination of integrators are needed. Only if this method selection does not work to your satisfaction, then I'd fiddle with giving specific integrators. What does NDSolve return for your specific system of ODEs?
– user21
Jun 19, 2013 at 7:11
• NDSolve returns an interpolating function(with the domain and range) which I can only plot, but it doesn't return the result as a function. Is there any command to get the function from an interpolating function. Jun 19, 2013 at 15:55
• @ruebenko I manually select the integrator a) for pedagogical purposes and b) to (try to) reproduce the results of someone else's integration scheme (e.g. matlab in this case). For my own stuff I let Mathematica handle it automatically. Also, sometimes I have to "port" to Fortran, where I may have LSODA at tops, sometimes less. I need to know how well the less sophisticated methods perform. Jun 19, 2013 at 16:50
• @user4402 This is a numerical solution to the ODE, not an analytic solution. If you want to create your own approximate representation, then you can extract the data from the InterpolatingFunction itself, or use it to generate the data at points that you need. Jun 19, 2013 at 17:04

Here is how to define a 5(4) pair of Dormand and Prince coefficients [DP80]. This is currently the method used by ode45 in MATLAB.

DOPRIamat = {
{1/5},
{3/40, 9/40},
{44/45, -56/15, 32/9},
{19372/6561, -25360/2187, 64448/6561, -212/729},
{9017/3168, -355/33, 46732/5247, 49/176, -5103/18656},
{35/384, 0, 500/1113, 125/192, -2187/6784, 11/84}};
DOPRIbvec = {35/384, 0, 500/1113, 125/192, -2187/6784, 11/84, 0};
DOPRIcvec = {1/5, 3/10, 4/5, 8/9, 1, 1};
DOPRIevec = {71/57600, 0, -71/16695, 71/1920, -17253/339200,
22/525, -1/40};
DOPRICoefficients[5, p_] :=
N[{DOPRIamat, DOPRIbvec, DOPRIcvec, DOPRIevec}, p];


Then:

NDSolve[system,
Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 5,
"Coefficients" -> DOPRICoefficients, "StiffnessTest" -> False}]


where

system


is the second order ODE specified in the usual Mathematica manner for ODE's.

• The simple NDSolve solves the ODE's i am using but it returns an interpolationg function. Whereas I want the function itself so that I can make further changes to it. How should I get the function from the interpolating function. Jun 19, 2013 at 16:07
• @user4402 -- Matlab's ODE45 returns a vector of values and the times at which those values occur. If you wish to duplicate this output, you can simply evaluate the interpolating function at the points you wish. For example, with fInt the interpolating function, Table[fInt[x],{x,0,10,0.1}] will return the values at the points 0, 0.1, 0.2, ... 10. Jun 19, 2013 at 17:21
• Thank you eric for the brilliant solution Apr 19, 2014 at 19:42

Check the documentation centre right here: LINK

And a simple example (modelling coupled springs with added nonlinear restoring forces (for large vibrations)):

eqn = {x''[t] ==
0.4 x[t] + -1/6 x[t]^3 - 1.808 (x[t] - y[t]) -
1/10 (x[t] - y[t])^3,
y''[t] == -1.808 (y[t] - x[t]) - 1/10 (y[t] - x[t])^3,
x[0] == -0.6, x'[0] == 1/2, y[0] == 3.001, y'[0] == 5.9};

sol = NDSolve[eqn, {x, y}, {t, 0, 200}]

ParametricPlot[Evaluate[{x[t], y[t]} /. sol], {t, 0, 200}]


ParametricPlot[Evaluate[{x[t], x'[t]} /. sol], {t, 0, 200}]


ParametricPlot[Evaluate[{y[t], y'[t]} /. sol], {t, 0, 200}]


• This is nice, but is it a coupled second-order system? I think that would be more helpful to the original poster. Jun 18, 2013 at 23:23
• Appreciate it ^^ Thank you ! Jun 20, 2013 at 10:11
• Thank you sector for your brilliant solution. Apr 19, 2014 at 19:44