--- a/libinterp/corefcn/dasrt.cc
+++ b/libinterp/corefcn/dasrt.cc
@@ -262,7 +262,6 @@
 of strings, inline functions, or function handles, the first element names
 the function @math{f} described above, and the second element names a
 function to compute the modified Jacobian
-
 @tex
 $$
 J = {\partial f \over \partial x}

--- a/libinterp/corefcn/dassl.cc
+++ b/libinterp/corefcn/dassl.cc
@@ -210,7 +210,6 @@
 of strings, inline functions, or function handles, the first element names
 the function @math{f} described above, and the second element names a
 function to compute the modified Jacobian
-
 @tex
 $$
 J = {\partial f \over \partial x}

--- a/libinterp/corefcn/gcd.cc
+++ b/libinterp/corefcn/gcd.cc
@@ -452,7 +452,6 @@
 
 Optional return arguments @var{v1}, @dots{}, contain integer vectors such
 that,
-
 @tex
 $g = v_1 a_1 + v_2 a_2 + \cdots$
 @end tex

--- a/libinterp/corefcn/givens.cc
+++ b/libinterp/corefcn/givens.cc
@@ -40,7 +40,6 @@
 @deftypefn  {} {@var{G} =} givens (@var{x}, @var{y})
 @deftypefnx {} {[@var{c}, @var{s}] =} givens (@var{x}, @var{y})
 Compute the Givens rotation matrix @var{G}.
-
 @tex
 The Givens matrix is a $2\times 2$ orthogonal matrix
 $$

--- a/libinterp/corefcn/gsvd.cc
+++ b/libinterp/corefcn/gsvd.cc
@@ -115,7 +115,6 @@
 
 The generalized singular value decomposition is defined by the following
 relations:
-
 @tex
 $$ A = U C X^\dagger $$
 $$ B = V S X^\dagger $$

--- a/libinterp/corefcn/ordqz.cc
+++ b/libinterp/corefcn/ordqz.cc
@@ -62,7 +62,6 @@
 Reorder the QZ@tie{}decomposition of a generalized eigenvalue problem.
 
 The generalized eigenvalue problem is defined as
-
 @tex
 $$A x = \lambda B x$$
 @end tex
@@ -96,7 +95,6 @@
 @var{ZR} such that the order of the eigenvalue on the diagonal of @var{AA} and
 @var{BB} is changed.  The resulting reordered matrices @var{AR} and @var{BR}
 fulfill:
-
 @tex
 $$ A_R = Q_R \cdot A \cdot Z_R, B_R = Q_R \cdot B \cdot Z_R $$
 @end tex

--- a/libinterp/corefcn/psi.cc
+++ b/libinterp/corefcn/psi.cc
@@ -49,7 +49,6 @@
 and so on.
 
 The digamma function is defined:
-
 @tex
 $$
 \Psi (z) = {d (log (\Gamma (z))) \over dx}

--- a/libinterp/corefcn/qz.cc
+++ b/libinterp/corefcn/qz.cc
@@ -66,7 +66,6 @@
 Compute the QZ@tie{}decomposition of a generalized eigenvalue problem.
 
 The generalized eigenvalue problem is defined as
-
 @tex
 $$A x = \lambda B x$$
 @end tex

--- a/libinterp/corefcn/schur.cc
+++ b/libinterp/corefcn/schur.cc
@@ -275,7 +275,6 @@
 complex, upper triangular Schur@tie{}form @var{T}.
 
 Note that the following relations hold:
-
 @tex
 $UR \cdot TR \cdot {UR}^T = U T U^{\dagger}$ and
 $U^{\dagger} U$ is the identity matrix I.

--- a/libinterp/corefcn/svd.cc
+++ b/libinterp/corefcn/svd.cc
@@ -83,7 +83,6 @@
 Compute the singular value decomposition of @var{A}.
 
 The singular value decomposition is defined by the relation
-
 @tex
 $$
  A = U S V^{\dagger}

--- a/scripts/general/cart2pol.m
+++ b/scripts/general/cart2pol.m
@@ -48,7 +48,6 @@
 ## @var{z}, if present, is unchanged by the transformation.
 ##
 ## The coordinate transformation is computed using:
-##
 ## @tex
 ## $$ \theta = \arctan \left ( {y \over x} \right ) $$
 ## $$ r = \sqrt{x^2 + y^2} $$
@@ -66,9 +65,6 @@
 ##
 ## @end ifnottex
 ##
-## @c FIXME: Remove this note in Octave 9.1 (two releases after 7.1).
-## Note: For @sc{matlab} compatibility, this function no longer returns a full
-## coordinate matrix when called with a single return argument.
 ## @seealso{pol2cart, cart2sph, sph2cart}
 ## @end deftypefn
 

--- a/scripts/general/cart2sph.m
+++ b/scripts/general/cart2sph.m
@@ -44,7 +44,6 @@
 ## @var{r} is the distance to the origin @w{(0, 0, 0)}.
 ##
 ## The coordinate transformation is computed using:
-##
 ## @tex
 ## $$ \theta = \arctan \left ({y \over x} \right ) $$
 ## $$ \phi = \arctan \left ( {z \over {\sqrt{x^2+y^2}}} \right ) $$
@@ -62,9 +61,6 @@
 ##
 ## @end ifnottex
 ##
-## @c FIXME: Remove this note in Octave 9.1 (two releases after 7.1).
-## Note: For @sc{matlab} compatibility, this function no longer returns a full
-## coordinate matrix when called with a single return argument.
 ## @seealso{sph2cart, cart2pol, pol2cart}
 ## @end deftypefn
 

--- a/scripts/general/divergence.m
+++ b/scripts/general/divergence.m
@@ -30,7 +30,6 @@
 ## @deftypefnx {} {@var{div} =} divergence (@var{fx}, @var{fy})
 ## Calculate divergence of a vector field given by the arrays @var{fx},
 ## @var{fy}, and @var{fz} or @var{fx}, @var{fy} respectively.
-##
 ## @tex
 ## $$
 ## div F(x,y,z) = \partial_x{F} + \partial_y{F} + \partial_z{F}

--- a/scripts/general/pol2cart.m
+++ b/scripts/general/pol2cart.m
@@ -47,7 +47,6 @@
 ## @var{z}, if present, is unchanged by the transformation.
 ##
 ## The coordinate transformation is computed using:
-##
 ## @tex
 ## $$ x = r \cos \theta $$
 ## $$ y = r \sin \theta $$
@@ -64,9 +63,7 @@
 ## @end example
 ##
 ## @end ifnottex
-## @c FIXME: Remove this note in Octave 9.1 (two releases after 7.1).
-## Note: For @sc{matlab} compatibility, this function no longer returns a full
-## coordinate matrix when called with a single return argument.
+##
 ## @seealso{cart2pol, sph2cart, cart2sph}
 ## @end deftypefn
 

--- a/scripts/general/rescale.m
+++ b/scripts/general/rescale.m
@@ -45,7 +45,6 @@
 ##
 ## Programming Notes:
 ## The applied formula is
-##
 ## @tex
 ## $$B = l + {A - inmin \over inmax - inmin} \cdot (u - l)$$
 ## @end tex

--- a/scripts/general/sph2cart.m
+++ b/scripts/general/sph2cart.m
@@ -45,7 +45,6 @@
 ## @var{r} is the distance to the origin @w{(0, 0, 0)}.
 ##
 ## The coordinate transformation is computed using:
-##
 ## @tex
 ## $$ x = r \cos \phi  \cos \theta $$
 ## $$ y = r \cos \phi  \sin \theta $$
@@ -62,9 +61,7 @@
 ## @end example
 ##
 ## @end ifnottex
-## @c FIXME: Remove this note in Octave 9.1 (two releases after 7.1).
-## Note: For @sc{matlab} compatibility, this function no longer returns a full
-## coordinate matrix when called with a single return argument.
+##
 ## @seealso{cart2sph, pol2cart, cart2pol}
 ## @end deftypefn
 

--- a/scripts/linear-algebra/planerot.m
+++ b/scripts/linear-algebra/planerot.m
@@ -27,7 +27,6 @@
 ## @deftypefn {} {[@var{G}, @var{y}] =} planerot (@var{x})
 ## Compute the Givens rotation matrix for the two-element column vector
 ## @var{x}.
-##
 ## @tex
 ## The Givens matrix is a $2\times 2$ orthogonal matrix
 ## $$

--- a/scripts/linear-algebra/vecnorm.m
+++ b/scripts/linear-algebra/vecnorm.m
@@ -31,7 +31,6 @@
 ## @var{dim}.
 ##
 ## The p-norm of a vector is defined as
-##
 ## @tex
 ## $$ {\Vert A \Vert}_p  = \left[ \sum_{i=1}^N {| A_i |}^p \right] ^ {1/p} $$
 ## @end tex

--- a/scripts/optimization/humps.m
+++ b/scripts/optimization/humps.m
@@ -29,7 +29,6 @@
 ## Evaluate a function with multiple minima, maxima, and zero crossings.
 ##
 ## The output @var{y} is the evaluation of the rational function:
-##
 ## @tex
 ## $$y = -{ {1200x^4 - 2880x^3 + 2036x^2 - 348x - 88} \over {200x^4 - 480x^3 + 406x^2 - 138x + 17} }$$
 ## @end tex

--- a/scripts/plot/draw/peaks.m
+++ b/scripts/plot/draw/peaks.m
@@ -32,7 +32,6 @@
 ## Plot a function with lots of local maxima and minima.
 ##
 ## The function has the form
-##
 ## @tex
 ## $$f(x,y) = 3 (1 - x) ^ 2 e ^ {\left(-x^2 - (y+1)^2\right)} - 10 \left({x \over 5} - x^3 - y^5\right) - {1 \over 3} e^{\left(-(x+1)^2 - y^2\right)}$$
 ## @end tex

--- a/scripts/polynomial/padecoef.m
+++ b/scripts/polynomial/padecoef.m
@@ -28,7 +28,6 @@
 ## @deftypefnx {} {[@var{num}, @var{den}] =} padecoef (@var{T}, @var{N})
 ## Compute the @var{N}th-order Pad@'e approximant of the continuous-time
 ## delay @var{T} in transfer function form.
-##
 ## @tex
 ## The Pad\'e approximant of $e^{-sT}$ is defined by the following equation
 ## $$ e^{-sT} \approx {P_n(s) \over Q_n(s)} $$

--- a/scripts/specfun/ellipke.m
+++ b/scripts/specfun/ellipke.m
@@ -42,7 +42,6 @@
 ## Mathematical Note:
 ##
 ## Elliptic integrals of the first kind are defined as
-##
 ## @tex
 ## $$
 ## {\rm K} (m) = \int_0^1 {dt \over \sqrt{(1 - t^2) (1 - m t^2)}}
@@ -63,7 +62,6 @@
 ## @end ifnottex
 ##
 ## Elliptic integrals of the second kind are defined as
-##
 ## @tex
 ## $$
 ## {\rm E} (m) = \int_0^1 {\sqrt{1 - m t^2} \over \sqrt{1 - t^2}} dt

--- a/scripts/specfun/expint.m
+++ b/scripts/specfun/expint.m
@@ -28,7 +28,6 @@
 ## Compute the exponential integral.
 ##
 ## The exponential integral is defined as:
-##
 ## @tex
 ## $$
 ## {\rm E_1} (x) = \int_x^\infty {e^{-t} \over t} dt

--- a/scripts/specfun/legendre.m
+++ b/scripts/specfun/legendre.m
@@ -43,7 +43,6 @@
 ## dimension more than @var{x}.
 ##
 ## The associated Legendre function of degree @var{n} and order @var{m}:
-##
 ## @tex
 ## $$
 ## P^m_n(x) = (-1)^m (1-x^2)^{m/2}{d^m\over {dx^m}}P_n (x)
@@ -63,7 +62,6 @@
 ##
 ## @noindent
 ## with Legendre polynomial of degree @var{n}:
-##
 ## @tex
 ## $$
 ## P(x) = {1\over{2^n n!}}\biggl({d^n\over{dx^n}}(x^2 - 1)^n\biggr)
@@ -101,7 +99,6 @@
 ## Legendre functions by the following:
 ##
 ## For Legendre functions of degree @var{n} and order 0:
-##
 ## @tex
 ## $$
 ## SP^0_n (x) = P^0_n (x)
@@ -120,7 +117,6 @@
 ## @end ifnottex
 ##
 ## For Legendre functions of degree n and order m:
-##
 ## @tex
 ## $$
 ## SP^m_n (x) = P^m_n (x)(-1)^m\biggl({2(n-m)!\over{(n+m)!}}\biggl)^{0.5}
@@ -144,7 +140,6 @@
 ## Legendre functions by the following:
 ##
 ## For Legendre functions of degree @var{n} and order @var{m}
-##
 ## @tex
 ## $$
 ## NP^m_n (x) = P^m_n (x)(-1)^m\biggl({(n+0.5)(n-m)!\over{(n+m)!}}\biggl)^{0.5}

--- a/scripts/statistics/kendall.m
+++ b/scripts/statistics/kendall.m
@@ -45,7 +45,6 @@
 ## is the correlation of the signs of all rank differences of
 ## @var{x} and @var{y}; i.e., if both @var{x} and @var{y} have distinct
 ## entries, then
-##
 ## @tex
 ## $$ \tau = {1 \over N(N-1)} \sum_{i,j} {\rm sign}(q_i-q_j) \, {\rm sign}(r_i-r_j) $$
 ## @end tex

--- a/scripts/statistics/moment.m
+++ b/scripts/statistics/moment.m
@@ -32,7 +32,6 @@
 ## Compute the @var{p}-th central moment of the vector @var{x}.
 ##
 ## The @var{p}-th central moment of @var{x} is defined as:
-##
 ## @tex
 ## $$
 ## {\sum_{i=1}^N (x_i - \bar{x})^p \over N}
@@ -87,7 +86,6 @@
 ##
 ## @item @qcode{"r"}
 ##   Raw Moment.  The moment about zero defined as
-##
 ## @tex
 ## $$
 ## {\rm moment} (x) = { \sum_{i=1}^N {x_i}^p \over N }