Gauss-Newton algorithm implementation

I am trying to implement a Gauss-Newton algorithm for non-linear optimisation. Unfortunately despite searching through the library and the internet I can't figure out quite how to use it in my case. As such I am hoping that I have come to the right place to ask for some assistance? Thank you to anyone willing to help me.

My problem consists of minimising a function $f(\vec x)$ where $\mathbb R^n \ni \vec x = (x_1, \dots ,x_n)$ and $f:\mathbb R^n\rightarrow \mathbb R$. Physically $f$ is a chi squared value $f=y-e(\vec x)$ where $y$ is the target value and $e(\vec x)$ is the test function evaluated at guess $\vec x$.

I begin by computing $f(\vec x+\delta \vec x)$ and $f(\vec x)$. Here $\vec x+\delta \vec x=x_1+\delta x_1,\dots , x_n+\delta x_n$. I then evaluate numerical derivatives using the following quotient formula,
\begin{equation}
\frac{\partial f(\vec x)}{\partial x_i} = \frac{f(x_1,\dots , x_i+\delta x_i,\dots ,x_n)-f(\vec x)}{\delta x_i}\qquad \forall i =1,\dots ,n
\end{equation}
I must then combine this into the Jacobian $\vec J$,
\begin{equation}
\vec J = \bigg(\frac{\partial f(\vec x)}{\partial x_1}\ \frac{\partial f(\vec x)}{\partial x_2} \cdots \frac{\partial f(\vec x)}{\partial x_n}\bigg)
\end{equation}
I can take the transpose of $\vec J$ denoted $\vec J^T$.

So far I can code the above (in C++). However The next part becomes confusing for me personally. In particular I do not know how to write the following matrix equation.
\begin{equation}
(\vec J \vec J^T)\delta \vec x = \vec Jf
\end{equation}
In fact I am not certain that this is the correct procedure?

Therefore, with numerical derivatives in place how can I obtain $\delta \vec x$ for the next iteration of the descent? Again many thanks for your time.

Solve the equation multiplying both terms of the equation for the inverse of the hessian matrix $\vec{J} \vec{J}^T$ on the left side:

$$
\delta \vec{x} = (\vec{J} \vec{J}^T)^{-1} \cdot \vec{J}f
$$

StackExchange Computer Science Q#54792, answer score: 2