Chebyshev net in 3D
The question asked is a special case of the following more general one: Let $$ \mathbb{W} = \{x^1 + \cdots + x^n = 0\} \subset \mathbb{R}^n. $$ Given an $n$-dimensional Riemannian manifold $M$ and a smooth map $f: \mathbb{W} \rightarrow M$, can $f$ be extended to a map $f: \mathbb{R}^n \rightarrow M$ such that $$ |\partial_1f|^2= \cdots = |\partial_nf|^2 = 1\ ? $$ Here, note that $\partial_if(x) \in T_{f(x)}M$, and the norm is taken with respect to the Riemannian metric on $M$.
The first remark, which I'm sure Anton and Robert already know, is that if $f$ happens to be an embedding (which is the case in the original question), then we can view it as a choice of all $n$ co-ordinate functions along a hypersurface in $M$ and the question is whether these functions can be extended to co-ordinate functions where the diagonal elements of the metric tensor written with respect to these co-ordinates are identically equal to $1$.
To determine what kind of PDE this system is, we can linearize it. Let $\dot{f}$ denote an infinitesimal variation of $f$. The linearized system is given by $$ 2\partial_if\cdot\partial_i\dot{f} = \dot{h}_i. $$ If $f$ is an immersion, then $\partial_1f(x), \ldots, \partial_nf(x)$ are always a basis of $\mathbb{R}^n$. Therefore, solving the linearized system for $\dot{f}$ is equivalent to solving for the functions $$ u_1 = \partial_1 f\cdot \dot{f}, \ldots, u_n = \partial_nf\cdot \dot{f}. $$ ``Differentiating by parts'', the linearized system can be rewritten as $$ 2\partial_i u_i - g^{jk}u_j\cdot\nabla_i\partial_kf = \dot{h}_i. $$ Note that the top order terms fully uncouple into ODE's in each of the co-ordinate directions. The coupling that occurs in the zero-th order terms prevents this system from being purely a system of ODE's. However, it is easy to check that this system is indeed a symmetric hyperbolic system, where initial smooth data posed on $\mathbb{W}$ can be extended to a unique smooth solution on all of $\mathbb{R}^n$. You can find more about this in a paper I wrote with DeTurck (Duke Math. J. Vol.51, No. 2 (1984), 243-260).
Using appropriate regularity estimates and the inverse function theorem with the appropriate Banach norm, this implies existence and uniqueness of a solution to the original nonlinear system on a tubular neighbhorhood of $\mathbb{W}$. I do not know whether there is a global solution; this requires a more careful analysis of this particular system.
ADDED: Actually, since the system is fully nonlinear, you cannot use the Banach space implicit function theorem to solve the system as given. You can, however, do one of two thing: Either "prolong" (which means roughly differentiate the system and add the partial derivatives of $f$ as unknown functions) the system into a quasilinear system or just use the Nash-Moser iteration argument. Either way, you still get what I describe above.
ADDED (in response to Robert's comment): Robert is right that if you specify $f$ along $\mathbb{W}$, the extension of $f$ is not uniquely determined, due to the full nonlinearity. Robert indicates that there are exactly 2 choices when $n = 3$. You have to do the linear algebra carefully along $\mathbb{W}$ to see this and count the number of possibilities in arbitrary dimension. I haven't done this yet.
ADDED: If along $\mathbb{W}$, you set $v = \frac{1}{n}(\partial_1 f+\cdots \partial_n f)$ and $u_i = \partial_if - v$, then the $u_1, \dots, u_n$ are tangential derivatives of $f$ along $\mathbb{W}$ and therefore given by the initial data. These vectors along $\mathbb{W}$ satisfy the equations $$ u_i\cdot v = \frac{1}{2}[1 - |u_i|^2 - \frac{1}{n}(1 - \sum_i |u_i|^2)] $$ $$ v\cdot v = \frac{1}{n}(1 - \sum_i |u_i|^2) $$ Clearly, a necessary condition for a solution is that $$ |u_1|^2 + \cdots + |u_n|^2 \le 1. $$ If strict inequality holds on $W$, I believe that there are always exactly two solutions on $W$. These translate into a corresponding inequality that the tangential derivatives of $f$ along $\mathbb{W}$ must satisfy. And if strict inequality holds, then there are exactly two ways to extend $f$ onto a tubular neighborhood of $\mathbb{W}$ such that $f$ defines a Chebyshev net.
Well, there's no problem if the map $f$ is real-analytic and close enough to the identity. Then it essentially reduces to the Cauchy-Kovalevskaya Theorem. For the smooth theory, you'll need something a bit more subtle, but it may be OK anyway. The characteristic variety (again, assuming that the initial conditions are sufficiently close to the identity) factors into three linear factors, which is something that Dennis DeTurck and Deane Yang know quite a bit about. It's possible that this problem is what they call `symmetric hyperbolic', in which case, the initial value problem will be OK. I'd check with one of them.
Oh: I should have explained that, because the characteristic variety is not smooth (being the union of three lines in general position in the projectivized cotangent space of each point of the domain, the initial value problem is not symmetric hyperbolic in the strict sense, so the simple hyperbolic existence theory for smooth initial data cannot be applied. That doesn't mean (as mentioned above) that a refinement of the symmetric hyperbolic theory won't work. The difference between this and the 2D case is that, in that case, the characteristic variety is 2 (real distinct, multiplicity one) points in projectivized cotangent space of each domain point, so, of course, it is smooth, and the standard hyperbolic existence theory applies.
By the way, you don't have uniqueness. When the obvious inequalities on the initial data are satisfied (so that the hyperbolic theory can be applied), there will be two (and only two) distinct solutions for each choice of initial data $f$ that you prescribe along $\mathbb{W}$.