Is there a numerical ODE solver facility in Maxima

(%i1) ? rk;
  -- Function: rk (<ODE>, <var>, <initial>, <domain>)
  -- Function: rk ([<ODE1>,...,<ODEm>], [<v1>,...,<vm>],
           [<init1>,...,<initm>], <domain>)
      The first form solves numerically one first-order ordinary
      differential equation, and the second form solves a system of m of
      those equations, using the 4th order Runge-Kutta method. <var>
      represents the dependent variable. <ODE> must be an expression
      that depends only on the independent and dependent variables and
      defines the derivative of the dependent variable with respect to
      the independent variable.
      The independent variable is specified with `domain', which must be
      a list of four elements as, for instance:
           [t, 0, 10, 0.1]
      the first element of the list identifies the independent variable,
      the second and third elements are the initial and final values for
      that variable, and the last element sets the increments that
      should be used within that interval.
      If <m> equations are going to be solved, there should be <m>
      dependent variables <v1>, <v2>, ..., <vm>. The initial values for
      those variables will be <init1>, <init2>, ..., <initm>.  There
      will still be just one independent variable defined by `domain',
      as in the previous case. <ODE1>, ..., <ODEm> are the expressions
      that define the derivatives of each dependent variable in terms of
      the independent variable. The only variables that may appear in
      those expressions are the independent variable and any of the
      dependent variables. It is important to give the derivatives
      <ODE1>, ..., <ODEm> in the list in exactly the same order used for
      the dependent variables; for instance, the third element in the
      list will be interpreted as the derivative of the third dependent
      variable.
      The program will try to integrate the equations from the initial
      value of the independent variable until its last value, using
      constant increments. If at some step one of the dependent
      variables takes an absolute value too large, the integration will
      be interrupted at that point. The result will be a list with as
      many elements as the number of iterations made. Each element in
      the results list is itself another list with <m>+1 elements: the
      value of the independent variable, followed by the values of the
      dependent variables corresponding to that point.
   There are also some inexact matches for `rk'.
   Try `?? rk' to see them.
(%o1) true

Best regards,
Martin