Let $\pi_n\colon\mathbb{S}^n\rightarrow\mathbb{R}P^n$ be the canonical projection, then $\pi_n$ is a smooth local diffeomorphism.
I am aware that I can directly conclude using that $\{\pm 1\}$ acts freely and properly discontinuously on the space $$\mathbb{R}P^n=\mathbb{S}^n/\{\pm 1\}$$ so that $\pi_n$ is a covering space. Besides, the smoothness of $\pi_n$ is easily read in charts.
However, I would like to read directly the desired property in charts using the inverse function theorem.
To do so, let $p\in\mathbb{S}^n$, one can assume without loss of generality that $p$ lies in the northern hemisphere so that: $$\phi\colon x\mapsto (x_1,\ldots,x_n)$$ is a chart of $\mathbb{S}^n$ around $p$ whose inverse is given by: $$\phi^{-1}\colon (x_1,\ldots,x_n)\mapsto\left(x_1,\ldots,x_n,\sqrt{1-{x_1}^2-\cdots-{x_n}^2}\right).$$ Furthermore, the following map is a chart of $\mathbb{R}P^n$ around $\pi_n(p)$: $$\psi\colon[x_1:\cdots:x_n:x_{n+1}]\mapsto \left(\frac{x_1}{x_{n+1}},\ldots,\frac{x_n}{x_{n+1}}\right).$$ Read in these above charts $\pi_n$ gives rise to: $$g:=\psi\circ\pi_n\circ\phi^{-1}\colon(x_1,\ldots,x_n)\mapsto\left(\frac{x_1}{\sqrt{1-{x_1}^2-\cdots-{x_n}^2}},\ldots,\frac{x_n}{\sqrt{1-{x_1}^2-\cdots-{x_n}^2}}\right).$$ Hence the smoothness of $\pi_n$ at $p$ and one has that:
$$\textrm{Jac}_p(g)_{i,j}=\frac{\delta_{i,j}}{\sqrt{1-{p_1}^2-\cdots-{p_n}^2}}-\frac{{p_i}^2}{(1-{p_1}^2-\cdots-{p_n}^2)^{3/2}}.$$
Believe me or not, I do not want to compute the nasty determinant of this matrix! Is there a clever way to see that the above matrix is invertible? Sadly, it does not appear to be diagonally strictly dominant. Any help will be greatly appreciated!
Let $r^2 = \sum x_i^2$. Then your map $g$ sends $(x_1,..., x_n)$ to $\frac{1}{\sqrt{1-r^2}}(x_1,..., x_n)$.
Now, the question is basically: given the (non-linear) system of equations $\frac{1}{\sqrt{1-r^2}} x_i = y_i$, can we solve for $x_i$?
Notice first that if we can somehow compute $r^2$ from this system solely in terms of the $y_i$ variables, then we can easily solve for $x_i$: $x_i = y_i \sqrt{1-r^2}$.
So lets do that. Square each of the equations, getting $\frac{x_i^2}{1-r^2} = y_i^2,$ and then add them all together. This gives $\frac{r^2}{1-r^2} = \sum y_i^2$. Clearing denominators gives $r^2 = \sum y_i^2 - r^2 \sum y_i^2$. Adding $r^2\sum y_i^2$ to both sides, factoring, and then dividing yields $$r^2 = \frac{\sum y_i^2}{1+\sum y_i^2}$$ and therefore, that $$x_i = y_i \sqrt{1-\frac{\sum y_i^2}{1+\sum y_i^2}}.$$
In other words, I am claiming that $g^{-1}$ maps $(y_1,...,y_n)$ to $\sqrt{1-\frac{\sum y_i^2}{1+\sum y_i^2}}(y_1,...., y_n)$.
Let's check it. We see that \begin{align*} (g^{-1}\circ g)(x_1,..., x_n) &= g^{-1}\left( \frac{x_1}{\sqrt{1-r^2}}, ..., \frac{x_n}{\sqrt{1-r^2}}\right)\\ &= \sqrt{1-\frac{ \frac{1}{1-r^2}\sum x_i^2}{1+ \frac{1}{1-r^2} \sum x_i^2}}\left(\frac{x_1}{\sqrt{1-r^2}}, ..., \frac{x_n}{\sqrt{1-r^2}}\right).\end{align*}
Focussing on that horrid coefficient and recalling that $r^2 = \sum x_i^2$, we see that it simplifies to $$\sqrt{1 - \frac{\frac{r^2}{1-r^2}}{1+ \frac{r^2}{1-r^2}}} = \sqrt{1 - \frac{\frac{r^2}{1-r^2}}{\frac{1}{1-r^2}}} = \sqrt{1-r^2}.$$
Thus, $(g^{-1}\circ g)(x_1,..., x_n) = (x_1,..., x_n)$ as claimed. You could check that $g\circ g^{-1}$ is the identity as well (but I leave that to you!), but note that you don't need to - From the chain rule, $Jac(g^{-1})Jac(g) = Id$, so your matrix must be invertible.