--- a/doc/interpreter/diffeq.txi	Mon May 18 16:12:23 2020 -0400
+++ b/doc/interpreter/diffeq.txi	Tue May 19 17:40:08 2020 +0900
@@ -120,14 +120,14 @@
 @node Matlab-compatible solvers
 @subsection Matlab-compatible solvers
 
-Octave also provides a set of solvers for initial value problems for Ordinary
-Differential Equations that have a @sc{matlab}-compatible interface.
+Octave also provides a set of solvers for initial value problems for ordinary
+differential equations (ODEs) that have a @sc{matlab}-compatible interface.
 The options for this class of methods are set using the functions.
 
 @itemize
-  @item @code{odeset}
+  @item @ref{XREFodeset,,odeset}
 
-  @item @code{odeget}
+  @item @ref{XREFodeget,,odeget}
 @end itemize
 
 Currently implemented solvers are:
@@ -136,37 +136,47 @@
   @item @nospell{Runge-Kutta} methods
 
   @itemize
-    @item @code{ode23} integrates a system of non-stiff ordinary differential
-    equations (ODEs) or index-1 differential-algebraic equations (DAEs).  It
-    uses the third-order @nospell{Bogacki-Shampine} method and adapts the local
-    step size in order to satisfy a user-specified tolerance.  The solver
-    requires three function evaluations per integration step.
+    @item @ref{XREFode45,,ode45} integrates a system of non-stiff ODEs or
+    index-1 differential-algebraic equations (DAEs) using the high-order,
+    variable-step @nospell{Dormand-Prince} method.  It requires six function
+    evaluations per integration step, but may take larger steps on smooth
+    problems than @code{ode23}: potentially offering improved efficiency at
+    smaller tolerances.
 
-    @item @code{ode45} integrates a system of non-stiff ODEs (or index-1 DAEs)
-    using the high-order, variable-step @nospell{Dormand-Prince} method.  It
-    requires six function evaluations per integration step, but may
-    take larger steps on smooth problems than @code{ode23}: potentially
-    offering improved efficiency at smaller tolerances.
+    @item @ref{XREFode23,,ode23} integrates a system of non-stiff ODEs or (or
+    index-1 DAEs).  It uses the third-order @nospell{Bogacki-Shampine} method
+    and adapts the local step size in order to satisfy a user-specified
+    tolerance.  The solver requires three function evaluations per integration
+    step.
+
+    @item @ref{XREFode23s,,ode23s} integrates a system of stiff ODEs (or
+    index-1 DAEs) using a modified second-order Rosenbrock method.
   @end itemize
 
   @item Linear multistep methods
 
   @itemize
-    @item @code{ode15s} integrates a system of stiff ODEs (or index-1 DAEs)
-    using a variable step, variable order method based on Backward Difference
-    Formulas (BDF).
+    @item @ref{XREFode15s,,ode15s} integrates a system of stiff ODEs (or
+    index-1 DAEs) using a variable step, variable order method based on
+    Backward Difference Formulas (BDF).
 
-    @item @code{ode15i} integrates a system of fully-implicit ODEs (or index-1
-    DAEs) using the same variable step, variable order method as @code{ode15s}.
-    The function @code{decic} can be used to compute consistent initial
-    conditions.
+    @item @ref{XREFode15i,,ode15i} integrates a system of fully-implicit ODEs
+    (or index-1 DAEs) using the same variable step, variable order method as
+    @code{ode15s}.  @ref{XREFdecic,,decic} can be used to compute consistent
+    initial conditions for @code{ode15i}.
   @end itemize
 @end itemize
 
+Detailed information on the solvers are given in @nospell{L. F. Shampine} and
+@nospell{M. W. Reichelt}, @cite{The MATLAB ODE Suite}, SIAM Journal on
+Scientific Computing, Vol. 18, 1997, pp. 1–22.
+
 @DOCSTRING(ode45)
 
 @DOCSTRING(ode23)
 
+@DOCSTRING(ode23s)
+
 @DOCSTRING(ode15s)
 
 @DOCSTRING(ode15i)

--- a/scripts/help/__unimplemented__.m	Mon May 18 16:12:23 2020 -0400
+++ b/scripts/help/__unimplemented__.m	Tue May 19 17:40:08 2020 +0900
@@ -73,8 +73,8 @@
       txt = ["Octave provides lsode and ode15i, ode15s, ode23, ode23s, ", ...
              "and ode45 for solving differential equations.  For more", ...
              "information try @code{help lsode}, @code{help ode45}.  ", ...
-             "Matlab-compatible ODE functions are provided by the odepkg ", ...
-             "package.  See @url{https://octave.sourceforge.io/odepkg/}."];
+             "Further specialized ODE solvers are provided by the odepkg ", ...
+             "package.  See @url{https://wiki.octave.org/Odepkg}."];
 
     case "polarplot"
       txt = ["polarplot is not implemented.  Consider using polar."];

--- a/scripts/ode/ode23s.m	Mon May 18 16:12:23 2020 -0400
+++ b/scripts/ode/ode23s.m	Tue May 19 17:40:08 2020 +0900
@@ -31,22 +31,21 @@
 ## @deftypefnx {} {@var{solution} =} ode23s (@dots{})
 ##
 ## Solve a set of stiff Ordinary Differential Equations (stiff ODEs) with a
-## Rosenbrock method of order (2,3).  For the definition of this method see
-## "The MATLAB ode suite", L.F. Shampine, mass.W. Reichelt, pp.6-7
+## Rosenbrock method of order (2,3).
 ##
 ## @var{fun} is a function handle, inline function, or string containing the
-## name of the function that defines the ODE: @code{M y' = f(t,y)}.  The function
-## must accept two inputs where the first is time @var{t} and the second is a
-## column vector of unknowns @var{y}. @var{M} is a costant mass matrix,
-## non-singular and possibly sparse. Set the field @var{Mass} in @var{odeopts}
-## using @var{odeset} to specify a mass matrix.
+## name of the function that defines the ODE: @code{M y' = f(t,y)}.  The
+## function must accept two inputs where the first is time @var{t} and the
+## second is a column vector of unknowns @var{y}. @var{M} is a constant mass
+## matrix, non-singular and possibly sparse. Set the field @var{Mass} in
+## @var{odeopts} using @var{odeset} to specify a mass matrix.
 ##
 ## @var{trange} specifies the time interval over which the ODE will be
-## evaluated.  Typically, it is a two-element vector specifying the initial and
-## final times (@code{[tinit, tfinal]}).  If there are more than two elements
-## then the solution will also be evaluated at these intermediate time
-## instances using an interpolation procedure of the same order as the one of
-## the solver.
+## evaluated.  Typically, it is a two-element vector specifying the initial
+## and final times (@code{[tinit, tfinal]}).  If there are more than two
+## elements then the solution will also be evaluated at these intermediate
+## time instances using an interpolation procedure of the same order as the
+## one of the solver.
 ##
 ## By default, @code{ode23s} uses an adaptive timestep with the
 ## @code{integrate_adaptive} algorithm.  The tolerance for the timestep
@@ -58,9 +57,9 @@
 ## the solution for the corresponding initial value in @var{init}.
 ##
 ## The optional fourth argument @var{ode_opt} specifies non-default options to
-## the ODE solver.  It is a structure generated by @code{odeset}.  @code{ode23s}
-## will ignore the following options: "BDF", "InitialSlope", "MassSingular",
-## "MStateDependence", "MvPattern", "MaxOrder", "NonNegative".
+## the ODE solver.  It is a structure generated by @code{odeset}.
+## @code{ode23s} will ignore the following options: "BDF", "InitialSlope",
+## "MassSingular", "MStateDependence", "MvPattern", "MaxOrder", "NonNegative".
 ##
 ## The function typically returns two outputs.  Variable @var{t} is a
 ## column vector and contains the times where the solution was found.  The