Approximating the implicit function Last week in our Calculus II course we learned about the implicit function theorem, to prove this theorem we used the inverse function theorem which we proved using using the Banach fixed-point theorem:
In the proof of the inverse function theorem we defined the following function and showed that it is a contraction: $$ \Phi_{y}:\,x\mapsto x+df(0)^{-1}(y-f(x)) $$
This provided us with a nice way to find the preimage of a function under a few assumptions, known as Newton iteration ($n=1$). Now concerning the implicit function theorem: I got the feeling that there might exist a similar way of approximating the implied function.
Let $ f\in C^1(\mathbb{R}^n;\mathbb{R}^l)$ and $p_0\in \mathbb{R}^n$ with $f(p_0)=0$. We choose coordinates $(x,y)\in\mathbb{R}^{n-l}\times\mathbb{R}^l\cong \mathbb{R}^n $, such that $det(\partial_yf(p_0))\neq0$. 
With the implicit function theorem there must exist a function $h\in C^1(U;\mathbb{R}^l)$, such that for all $x\in U: f(x, h(x))=0$. Where $U$ is a neighboorhood of $x_0\in \mathbb{R}^{n-l}$ with $p_0=(x_0,y_0)$.
My question is: Suppose we have a point $p_0$ and a function $f$ as mentioned above, is there a way of (iteratively) approximate $h(x)$ for any $x\in U$? Is there even a contraction which approximates $h$ itself (i.e. some function that operates in the $C^1$ space with supremum norm)?
Note: Since we can calulate $dh(x)=-(\partial_yf(x,h(x)))^{-1}\cdot\partial_xf(x,h(x))$, we can find a linear approximation. But I'm looking for an approximation which we can tune in terms of precision.
 A: Suppose we have some function like in the question, but let's keep it general: So $f(p_0) = s$ with some regular point $p_0=(x_0,y_0)$ in the chosen coordinates.
We define $F(x,y):=\left(\begin{array}{c}
x\\f(x,y)\end{array}\right)$ where $y=h(x)$. We know $dF(p_0)=\left(\begin{array}{cc}
\mbox{id} & 0\\
\frac{\partial f}{\partial x}(p_0) & \frac{\partial f}{\partial y}(p_0)
\end{array}\right)$ and thus we can find its inverse: $\left(\begin{array}{cc}
\mbox{id} & 0\\* & (\frac{\partial f}{\partial y}(0))^{-1}
\end{array}\right)$.
Now we can use the Newton-method to find the preimage of $\left(\begin{array}{c}
x\\s\end{array}\right)$ under F, by defining: 
$$ \Phi_{\left(\begin{array}{c}
x\\s\end{array}\right)}(x,y):=\left(\begin{array}{c}
x\\y\end{array}\right)+(dF(p_0))^{-1}(\left(\begin{array}{c}
x\\s\end{array}\right)-F(x,y))=\left(\begin{array}{c}
x\\y\end{array}\right)+\left(\begin{array}{cc}
\mbox{id} & 0\\* & (\frac{\partial f}{\partial y}(0))^{-1}
\end{array}\right)\cdot\left(\begin{array}{c}
0\\s-f(x,y)\end{array}\right)$$
Since we're only interested in the "second" component, we let
$$\phi_{\left(\begin{array}{c}
x\\s\end{array}\right)}:=y+(\frac{\partial f}{\partial y}(0))^{-1}(s-f(x,y)).$$
Now we have an approximating scheme for the implied function $h$:
$$h_0(x)=y_0$$
$$h_{i+1}(x)=\phi_{\left(\begin{array}{c}
x\\s\end{array}\right)}(h_i(x))=h_i(x)+(\frac{\partial f}{\partial y}(0))^{-1}(s-f(x,h_i(x)))$$

For example: If $f(x,y)=x^2+y^2$, let's try to find the solutions of $f(x,y)=1$, knowing that $p_0:=(0,1)$ is one:
$df(x,y)=(2x,\,\,\,2y)$ and $(\frac{\partial f}{\partial y}(0))^{-1}=\frac{1}{2}$
$$h_0(x)=1$$
$$h_1(x)=1+\frac{1}{2}(1-x^2-1)=\frac{1}{2}(2-x^2)$$
$$h_2(x)=\frac{1}{2}(2-x^2)+\frac{1}{2}(1-x^2-\frac{1}{4}(2-x^2)^2)=\frac{1}{8}(8-4x^2-x^4)$$
$$\dots$$
