diff --git a/libRALFit/doc/C/howtouse.rst b/libRALFit/doc/C/howtouse.rst
index 4eccaa95..7412312c 100644
--- a/libRALFit/doc/C/howtouse.rst
+++ b/libRALFit/doc/C/howtouse.rst
@@ -191,13 +191,13 @@ It must have the following signature:
 
    :param x: |eval_J_X|
 
-   :param J: |eval_J_r|
+   :param J: |eval_J_r| ``J(i*n+j)`` must be set to hold :math:`\nabla_{x_{j+1}} r_{i+1}( {\bm x} _{k}^{})`.
 
    :param status: |eval_J_status|
 
 
-For evaluating the function :math:`Hf = \sum_{i=1}^m r_i( {\bm x} )  {\bm W} \nabla^2 r_i( {\bm x} )`
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+For evaluating the function :math:`Hf = \sum_{i=1}^m r_i( {\bm x} )  \nabla^2 r_i( {\bm x} )`
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 A subroutine must be supplied to calculate
 :math:`Hf = \sum_{i=1}^m ( {\bm r} )_i \nabla^2 r_i( {\bm x} )` for
diff --git a/libRALFit/doc/Fortran/howtouse.rst b/libRALFit/doc/Fortran/howtouse.rst
index e496882f..bc3ac7ab 100644
--- a/libRALFit/doc/Fortran/howtouse.rst
+++ b/libRALFit/doc/Fortran/howtouse.rst
@@ -207,12 +207,12 @@ A subroutine must be supplied to calculate
 
    :p real X(n) [in]: |eval_J_X|
 
-   :p real J(m*n) [out]: |eval_J_r|
+   :p real J(m*n) [out]: |eval_J_r| ``J((i-1)*n+j)`` must be set to hold :math:`\nabla_{x_j} r_i(  {\bm x} _{k}^{})`  if the option ``Fortran_Jacobian=.true.``.
 
    :p params_base_type params [in]: |eval_J_params|
 
 
-For evaluating the function :math:`Hf = \sum_{i=1}^m r_i( {\bm x} )  {\bm W} \nabla^2 r_i( {\bm x} )`
+For evaluating the function :math:`Hf = \sum_{i=1}^m r_i( {\bm x} )  \nabla^2 r_i( {\bm x} )`
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 A subroutine must be supplied to calculate
@@ -253,11 +253,11 @@ subroutine must implement the following interface:
 
    :p params_base_type params [in]: |eval_Hf_params|
 
-For evaluating the function :math:`P({\bm x},{\bm y}) := ( H_1({\bm x}){\bm y} \dots  H_m({\bm x}){\bm y})`
+For evaluating the function :math:`P({\bm x},{\bm y}) := (  \nabla^2 r_1({\bm x}){\bm y} \dots   \nabla^2 r_m({\bm x}){\bm y})`
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 A subroutine may be supplied to calculate
-:math:`P({\bm x},{\bm y}) := ( H_1({\bm x}){\bm y} \dots  H_m({\bm x}){\bm y})` for
+:math:`P({\bm x},{\bm y}) := ( \nabla^2 r_1({\bm x}){\bm y} \dots  \nabla^2 r_m({\bm x}){\bm y})` for
 given vectors :math:`{\bm x}, {\bm y} \in \mathbb{R}^n`. The
 subroutine must implement the following interface:
 
diff --git a/libRALFit/doc/common/subroutines.rst b/libRALFit/doc/common/subroutines.rst
index b05066ce..3628d96b 100644
--- a/libRALFit/doc/common/subroutines.rst
+++ b/libRALFit/doc/common/subroutines.rst
@@ -18,7 +18,7 @@
 
 .. |weights| replace:: If present, this holds the square-roots of the diagonal entries of the weighting matrix, :math:`{\bm W}`. If absent, then the norm in the least squares problem is taken to be the 2-norm, that is, :math:`{\bm W} = I`.
 
-.. |eval_HP_desc| replace:: If present, this is a routine that, given vectors :math:`{\bm x}, {\bm y} \in \mathbb{R}^m`, returns the matrix :math:`P({\bm x},{\bm y}) := ( H_1({\bm x}){\bm y} \dots  H_m({\bm x}){\bm y})`. Further details of the format required are given in |eval_HP| in :ref:`user-routines`.  This is only referenced if ``model = 4`` in |nlls_options|.
+.. |eval_HP_desc| replace:: If present, this is a routine that, given vectors :math:`{\bm x}, {\bm y} \in \mathbb{R}^m`, returns the matrix :math:`P({\bm x},{\bm y}) := ( \nabla^2 r_1({\bm x}){\bm y} \dots  \nabla^2 r_m({\bm x}){\bm y})`. Further details of the format required are given in |eval_HP| in :ref:`user-routines`.  This is only referenced if ``model = 4`` in |nlls_options|.
 		       
 .. |iterate_X| replace:: on the first call it must hold the initial guess for :math:`\bm x`. On return it holds the value of :math:`\bm x` at the current iterate,  and must be passed unaltered to any subsequent call to |nlls_iterate|.
 
@@ -44,7 +44,7 @@
 
 .. |eval_J_X| replace::  holds the current point :math:`{\bm x}_{k}^{}` at which to evaluate :math:`{\bm J} (  {\bm x} _{k}^{})`.
 
-.. |eval_J_r| replace:: must be set by the routine to hold the Jacobian of the residual function evaluated at the current point :math:`{\bm x}_{k}^{}`, :math:`{\bm r} (  {\bm x} _{k}^{})`. ``J(i*m+j)`` must be set to hold :math:`\nabla_{x_j} r_i(  {\bm x} _{k}^{})`.
+.. |eval_J_r| replace:: must be set by the routine to hold the Jacobian of the residual function evaluated at the current point :math:`{\bm x}_{k}^{}`, :math:`{\bm r} (  {\bm x} _{k}^{})`. 
 
 .. |eval_J_status| replace:: is initialised to ``0`` before the routine is called. If it is set to a non-zero value by the routine, then |nlls_solve| / |nlls_iterate| will exit with error.
 
@@ -72,6 +72,6 @@
 
 .. |eval_HP_y| replace:: holds :math:`{\bm y}`, the vector which multiplies each Hessian.
 
-.. |eval_HP_HP| replace:: must be set by the routine to hold the matrix :math:`P({\bm x},{\bm y}) := ( H_1({\bm x}){\bm y} \dots  H_m({\bm x}){\bm y})`, held by columns as a vector.
+.. |eval_HP_HP| replace:: must be set by the routine to hold the matrix :math:`P({\bm x},{\bm y}) := ( \nabla^2 r_1({\bm x}){\bm y} \dots   \nabla^2 r_m({\bm x}){\bm y})`, held by columns as a vector.
 
 .. |eval_HP_status| replace:: is initialised to ``0`` before the routine is called. If it is set to a non-zero value by the routine, then |nlls_solve| / |nlls_iterate| will exit with error.