--- a/NEWS	Mon Jun 07 20:43:13 2021 +0200
+++ b/NEWS	Tue Jun 08 14:27:27 2021 +0900
@@ -61,6 +61,7 @@
 
 ### Documentation
 
+- Improve Differential Equations chapter and example for lsode (bug #57508).
 - Clarify usage of "Depends" keyword in package `DESCRIPTION` file (bug #60053).
 - Add note that wildcard patterns for `save` are glob patterns (bug #60157).
 - Change example for Delaunay triangulation to match the generating code (bug #60216).

--- a/doc/interpreter/diffeq.txi	Mon Jun 07 20:43:13 2021 +0200
+++ b/doc/interpreter/diffeq.txi	Tue Jun 08 14:27:27 2021 +0900
@@ -19,13 +19,13 @@
 @node Differential Equations
 @chapter Differential Equations
 
-Octave has built-in functions for solving ordinary differential equations,
-and differential-algebraic equations.
-All solvers are based on reliable ODE routines written in Fortran.
+Octave has built-in functions for solving ordinary differential equations
+(ODEs), and differential-algebraic equations (DAEs).
 
 @menu
 * Ordinary Differential Equations::
 * Differential-Algebraic Equations::
+* Matlab-compatible solvers::
 @end menu
 
 @cindex differential equations
@@ -104,21 +104,63 @@
 @end group
 @end example
 
-See @nospell{Alan C. Hindmarsh},
-@cite{ODEPACK, A Systematized Collection of ODE Solvers},
-in Scientific Computing, @nospell{R. S. Stepleman}, editor, (1983)
-for more information about the inner workings of @code{lsode}.
-
 An m-file for the differential equation used above is included with the
 Octave distribution in the examples directory under the name
 @file{oregonator.m}.
 
-@menu
-* Matlab-compatible solvers::
-@end menu
+
+@node Differential-Algebraic Equations
+@section Differential-Algebraic Equations
+
+The function @code{daspk} can be used to solve DAEs of the form
+@tex
+$$
+ 0 = f (\dot{x}, x, t), \qquad x(t=0) = x_0, \dot{x}(t=0) = \dot{x}_0
+$$
+@end tex
+@ifnottex
+
+@example
+0 = f (x-dot, x, t),    x(t=0) = x_0, x-dot(t=0) = x-dot_0
+@end example
+
+@end ifnottex
+
+@noindent
+where
+@tex
+$\dot{x} = {dx \over dt}$
+@end tex
+@ifnottex
+@math{x-dot}
+@end ifnottex
+is the derivative of @math{x}.  The equation is solved using
+@nospell{Petzold's} DAE solver @sc{daspk}.
+
+@DOCSTRING(daspk)
+
+@DOCSTRING(daspk_options)
+
+Octave also includes @sc{dassl}, an earlier version of @sc{daspk},
+and @sc{dasrt}, which can be used to solve DAEs with constraints
+(stopping conditions).
+
+@DOCSTRING(dassl)
+
+@DOCSTRING(dassl_options)
+
+@DOCSTRING(dasrt)
+
+@DOCSTRING(dasrt_options)
+
+See @nospell{K. E. Brenan}, et al., @cite{Numerical Solution of Initial-Value
+Problems in Differential-Algebraic Equations}, North-Holland (1989),
+DOI: @url{https://doi.org/10.1137/1.9781611971224},
+for more information about the implementation of @sc{dassl}.
+
 
 @node Matlab-compatible solvers
-@subsection Matlab-compatible solvers
+@section Matlab-compatible solvers
 
 Octave also provides a set of solvers for initial value problems for ordinary
 differential equations (ODEs) that have a @sc{matlab}-compatible interface.
@@ -169,7 +211,8 @@
 
 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.
+Scientific Computing, Vol. 18, 1997, pp. 1–22,
+DOI: @url{https://doi.org/10.1137/S1064827594276424}.
 
 @DOCSTRING(ode45)
 
@@ -188,51 +231,3 @@
 @DOCSTRING(odeget)
 
 @DOCSTRING(odeplot)
-
-@node Differential-Algebraic Equations
-@section Differential-Algebraic Equations
-
-The function @code{daspk} can be used to solve DAEs of the form
-@tex
-$$
- 0 = f (\dot{x}, x, t), \qquad x(t=0) = x_0, \dot{x}(t=0) = \dot{x}_0
-$$
-@end tex
-@ifnottex
-
-@example
-0 = f (x-dot, x, t),    x(t=0) = x_0, x-dot(t=0) = x-dot_0
-@end example
-
-@end ifnottex
-
-@noindent
-where
-@tex
-$\dot{x} = {dx \over dt}$
-@end tex
-@ifnottex
-@math{x-dot}
-@end ifnottex
-is the derivative of @math{x}.  The equation is solved using
-@nospell{Petzold's} DAE solver @sc{daspk}.
-
-@DOCSTRING(daspk)
-
-@DOCSTRING(daspk_options)
-
-Octave also includes @sc{dassl}, an earlier version of @sc{daspk},
-and @sc{dasrt}, which can be used to solve DAEs with constraints
-(stopping conditions).
-
-@DOCSTRING(dassl)
-
-@DOCSTRING(dassl_options)
-
-@DOCSTRING(dasrt)
-
-@DOCSTRING(dasrt_options)
-
-See @nospell{K. E. Brenan}, et al., @cite{Numerical Solution of Initial-Value
-Problems in Differential-Algebraic Equations}, North-Holland (1989) for
-more information about the implementation of @sc{dassl}.

--- a/libinterp/corefcn/lsode.cc	Mon Jun 07 20:43:13 2021 +0200
+++ b/libinterp/corefcn/lsode.cc	Tue Jun 08 14:27:27 2021 +0900
@@ -258,6 +258,22 @@
 
 You can use the function @code{lsode_options} to set optional
 parameters for @code{lsode}.
+
+See @nospell{Alan C. Hindmarsh},
+@cite{ODEPACK, A Systematized Collection of ODE Solvers},
+in Scientific Computing, @nospell{R. S. Stepleman}, editor, (1983)
+or @url{https://computing.llnl.gov/projects/odepack}
+for more information about the inner workings of @code{lsode}.
+
+Example: Solve the @nospell{Van der Pol} equation
+
+@example
+@group
+fvdp = @@(@var{y},@var{t}) [@var{y}(2); (1 - @var{y}(1)^2) * @var{y}(2) - @var{y}(1)];
+@var{t} = linspace (0, 20, 100);
+@var{y} = lsode (fvdp, [2; 0], @var{t});
+@end group
+@end example
 @seealso{daspk, dassl, dasrt}
 @end deftypefn */)
 {