I was reading a paper, where the author gives a method for solving a differential equations system using finite difference method. I am trying to simulate this result. The problem I am facing is that the author has mentioned mean value of Jacobian, which I am not sure what it means. I am attaching the screen shot of the theorem. The differential system is given as: $$\dot{x}=f(x,u,\bar{r}).$$ But you can assume that $\bar{r}$ and $u$ are constants.
Note: $x$ is a vector and $f$ is a vector valued function.
The use of term "mean value" is a reference to the underlying mean value theorem.
Specifically, let $f : \mathbb{R}^n \rightarrow \mathbb{R}$ be differentiable, and let $x, y \in \mathbb{R}^n$. Let $\phi : \mathbb{R} \rightarrow \mathbb{R}^n$ be given by \begin{equation} \phi(t) = f( (1-t) x + t y). \end{equation} Let $Df(x) \in \mathbb{R}^{n \times n}$ denote the Jacobian of $f$ at the point $x$. Then by the main theorem of calculus and the chain rule of differentiation \begin{multline} f(x) - f(y) = \phi(1) - \phi(0) = \int_0^1 \phi'(s) dx \\ = \int_0^1 \left[ Df((1-s)x + s y)\right]\cdot (x-y) dx = \left[ \int_0^1 Df((1-s)x + s y)ds\right] \cdot (x-y) = \overline{Df} \cdot (x - y), \end{multline} where we have the identification \begin{equation} \overline{Df} = \int_0^1 Df((1-s)x + s y)ds. \end{equation} One could argue that this matrix represents an average of the Jacobian over the line segment from $x$ to $y$. If $n = 1$ then the mean value theorem of integration implies that \begin{equation} \overline{Df} = Df((1-\xi) x + \xi y) \end{equation} for at least one $\xi \in (0,1)$, but this result does not carry to $n>1$ and it does not help us anyway.
From the point of view of practical calculations this expression is nearly useless and the author is only using it to prove that there exists a time step $\tau_j$ such that the $\|f^{(j+1)}\| < \|f^{j}\|$. It seems to me that a more practical approach would be to consider $f^{(j+1)}$ as a function of the time step, i.e. \begin{equation} f^{j+1} = f^{j+1}(\tau) = f(x^j + \tau \cdot f (x^j,u^j,\overline{r})) \end{equation} and then seek a value of $\tau^j \in \mathbb{R}$ which minimizes $\|f^{j+1}(\tau)\|^2$.
So define a function $\phi_j(\tau) = \|f^{j+1}(\tau)\|^2$ and solve the equation $\phi_j'(\tau) = 0$. I imagine that it will be hard to apply Newton's method as the expression for $\phi_j'$ is already complicated and $\phi''(\tau)$ will be worse. I think, you will have more luck with the secant method.
Keep in mind that the paper that you are reading, i.e. http://www.jstor.org/stable/pdf/1911749.pdf is from 1965. It is likely in the extreme that better methods have been developed in the interrim. That being said, there is nothing wrong in implementing old methods. I think that the the only way to truly understand the limitations of an algorithm is to implement it and experiment with it.