How calculate programs/packages like Matlab, Python/scipy, ...the inverse jacobian for multidimensional Newton's method? x_(n+1)=xn−(J(x_n)^(−1)∗f(x_n))

Most numerical packages give you the option of either computing the Jacobian yourself and passing it to the solver, or of numerically approximating it with a finite difference scheme. I imagine that in general while performing Newton's method or other methods expressed in terms of an inverse Jacobian, these packages do not actually compute the inverse for reasons of stability. Instead, they solve the linear system $J(x_n)x_{n+1}=J(x_n)x_n-f(x_n)$ for $x_{n+1}$ at each time step.
The Jacobian is not always invertible; in order to use Newton's method the Jacobian must be invertible, though. You can see this in one dimension, with a function $\mathbb{R}\to \mathbb{R}$ with a critical point which is not an extremum.