Lipschitz Constant for Equation System

Question

I'm would like to compute the region of convergence from Newton-Rapson Method for a multivariate system of equations by a theorem that states:

Let $\vec{F}: \mathbb{R}^N \rightarrow \mathbb{R}^N$ be continuously differentiable in an open convex set $D \subset \mathbb{R}^N$. Assume

1. there exists $x^*$ with $r,\beta > 0$ such that $N(x*,r) \subset D$. ($N$ ist the disk around $x^*$ with radius $r$)
2. $\vec{F}(x^*)=0$
3. $\Vert J^{-1} \Vert \leq \beta$
4. and J is Lipschitz continuous with Lipschitz constant $\gamma$

Then there exists $\epsilon>0$ such that for all $x_0 \in N(x*,r)$ the sequence $x_1, x_2 ,...$ converges to $x^*$ by applying the Newton-Rapson Method and obeys:

1. $\Vert x_{k+1} - x^* \Vert \leq \beta \gamma \Vert x_{k} - x^* \Vert^2$

I found an example within a book. Within this book, there is a Jacobian J $: \begin{bmatrix} 2 x_1 & 2 x_2\\ -1 & 2 x_2\\ \end{bmatrix}$ given to the underlying equation system F $: \begin{bmatrix} x_1^2 x_2^2 - 9\\ -x_1 x_2^2 - 3 \\ \end{bmatrix} = 0$ and a solution $x^*= [2, \sqrt5]$. The Lipschitz constant is said to be $\gamma=4$.

How can one compute the Lipschitz constant in general for multivariate system of equations and how do they determine $\gamma=4$?