I have a function $f \colon R^3 \to R$. I want to find $x$, $y$, $z$ such that $f(x,y,z)=0$.
I'm using the method from here: http://en.wikipedia.org/wiki/Quasi-Newton_method#Search_for_zeroes. Specifically, I need to find the Jacobian matrix of $f$. This is where I start getting confused, because the wikipedia page on Jacobian matrices says that it operates on a vector-valued function, while my function is scalar-valued. Maybe it's a vector of one element? Assuming this, I could go ahead and create a 1x3 Jacobian matrix, but the very next step is to invert the matrix, which is not defined for non-square matrices.
What am I misunderstanding here?
To get isolated points as solutions you need as many equations as there are variables. Each single equation in $n$ variables defines a hypersurface of co-dimension 1 o dimension $n-1$ (in ideal circumstances). The regular, transversal intersection of $k$ such hypersurfaces leads to a manifold of co-dimension $k$ or dimension $n-k$, an isolated point has dimension $0$ and is thus the regular intersection of $n$ hypersurfaces.
This should all be known in the case the equations are linear and thus the hypersurfaces and intersections affine subspaces.
You have one equation for 3 variables. In principle you could fix any two variables at an arbitrary value and solve for the third. In actuality, this rarely works, since there will be large areas that are not covered by the projection of the solution surface onto coordinate planes.
Think for example about $f(x,y,z)=x^2+y^2+z^2-1$ which only has for $x_0^2+y_0^2\le 1$ solutions $z$ for $f(x_0,y_0,z)=0$.
You could still follow the negative gradient of $\ln|f(x,y,z)|$ to fall into the singularity at the solution surface, i.e., solve the differential equations \begin{align} \dot x&=-\frac{\partial_xf(x,y,z)}{f(x,y,z)}\\ \dot y&=-\frac{\partial_yf(x,y,z)}{f(x,y,z)}\\ \dot z&=-\frac{\partial_zf(x,y,z)}{f(x,y,z)}\\ \end{align}