Diophantine equation to characterize natural numbers Let's consider
$$\Bbb N=\{0,1,2,3,\ldots\},$$
and, for each $k\in\{1,2,3,\ldots,\}$, let
$$o_k=2k-1$$
be the sequence of odd natural numbers.
Given that for each $m\in\Bbb N$, if $a$ is odd, the number
$$m(m+a)$$
is even, it is easy to see that
$$\begin{align*}
(x+y)^2 + y + 3x &= x^2 + 2xy + y^2 + y + 3x\\
&= x(x+3) + y(y+1) + 2xy
\end{align*}$$
is even.
Assuming that for arbitrary $x_1,\ldots, x_k\in\Bbb N$, the number
$$(x_1+\cdots + x_k)^2+x_k+3x_{k-1}+\cdots +o_kx_1$$
is even, we have that
$$\begin{align*}
& (x_1+\cdots + x_{k+1})^2+x_{k+1}+3x_k+\cdots +o_kx_1\\
&\qquad = (x_1+\cdots + x_{k})^2 +2(x_1+\cdots + x_{k})x_{k+1}+ x_{k+1}^2 + x_{k+1} + (x_k+\cdots + o_kx_{1}) +2(x_1+\cdots + x_{k}) \\
&\qquad = (x_1+\cdots + x_{k})^2 + (x_k+\cdots + o_kx_{1}) + x_{k+1}(x_{k+1}+1) + 2(x_1+\cdots + x_{k})(x_{k+1}+1)
\end{align*}$$
is even.
Therefore by induction, for each $k\geq 2$ we have that
$$\frac1{2}\left[(x_1+\cdots + x_k)^2+x_k+3x_{k-1}+\cdots +o_kx_1\right]$$
is indeed a natural number.
But, does this describes all natural numbers?
I mean, for a given $n\in\Bbb N$, does the equation
$$(x_1+\cdots + x_k)^2+x_k+3x_{k-1}+\cdots +o_kx_1=2n$$
have exactly one solution $(x_1,\ldots,x_k)\in\Bbb N^k$?
 A: The "Cantor pairing function" is $f(x,y)=(1/2)[(x+y)^2+3x+y]$ and is known to give a one-to-one correspondence between $N \times N$ and $N$. This is the function of the OP for $k=2$.
It is a conjecture that this $f(x,y)$ and it's flip $f(y,x)$ are the only polynomials with real coefficients that give such a one-to-one correspondence.
There is a result, the Fueter-Polya Theorem, that says these are the only quadratic polynomials with real coefficients that, when restricted to $N \times N$, map it one-to-one onto $N$. The proof is quite involved, see for example Craig Smorynski's book "Logical Number Theory I", chapter 1, sections 3,4,5. 
Of course it is relatively easy to verify these pairing functions do what one wants, the difficulty is showing they are the only ones that work, among polynomials.
As I noted in a comment, it seems your attempt to extend this doesn't quite work. That is, you seem to be extending the Cantor pairing function, in the specific case of $N \times N \times N$, to the function
$$f(x,y,z)=(1/2)[(x+y+z)^2+z+3y+5x].$$
However with this $f$ we have $f(1,0,0)=f(0,0,2)=6/2=3$. 
And also $f(0,0,3)=f(1,1,0)=6,$ and probably a lot of other collisions.
It would indeed be interesting, in my opinion, if the function could be adjusted even to the three variable case. Of course one can soup one up by defining $g(x,y,z)=f(f(x,y),z)$ (or something like that) but that is the same as one of the usual definitions of ordered triples via a previous definition of ordered pairs.
ADDED: There is an "order of magnitude" problem with trying to represent $n$-tuples distinctly by a quadratic function. Take for example the case of representing triples. Consider all triples $x,y,z$ with $x,y,z \le N$. There are $(N+1)^3$ of these triples. No quadratic function can work to distinctly represent these. For example if $x,y,z \le N$ then the function $f(x,y,z)=(1/2)[(x+y+z)^2+z+3y+5x]$ is bounded above by $(1/2)[(3N)^2+9N]$ and as $N \to \infty$ cannot represent distinctly $(N+1)^3$ triples. This order of magnitude problem gets worse for $k$-tuples as $k=4,5,6...$.
I think to represent $k$- tuples one would need the polynomial to be of degree at least $k$ because of order of magnitude issues.
A Stab at triples (and beyond): Let $C(m,k)$ denote the binomial coefficient, where we understand that if $m<k$ then $C(m,k)=0.$ Then the following (I think) maps the triples in $N^3$ one-to-one onto $N$:
$$f(x,y,z)=C(x+y+z+2,3)+C(x+y+1,2)+C(x,1).$$
Note that this is in a way a natural extension of the Cantor pairing, since one form of that pairing is $C(x+y+1,2)+C(x,1).$
For 4-tuples the idea would be
$$f(w,x,y,z)=C(w+x+y+z+3,4)+C(w+x+y+2,3)+C(w+x+1,2)+C(w,1),$$
and generally for $k$-tuples the pattern would start with a coefficient $C(s+k-1,k)$ where $s$ is the sum of the tuple, then $C(t+k-2,k-1)$ where $t$ is the sum of the first $k-1$ entries of the $k$-tuple, and so on.
I haven't (yet) come up with a proof that these maps indeed take the $k$-tuples one-to-one onto $N$, but have checked that the $k=3$ case works for an initial segment of triples, by going through the triangular parts $x+y+z=const$ one at a time, enumerating each in turn.
A: Since $f(x,y)=\frac12[(x+y)^2+3x+y]$ is a bijection $\Bbb N\times\Bbb N\leftrightarrow\Bbb N$, we get that $g(u,v,w)=f(f(u,v),w)$ is a bijection $\Bbb N^3\leftrightarrow\Bbb N$. However, $g$ is of degree $4$ for $3$ variables.
Still, notice that for $k=2^r$ we can construct a bijection $h_r:\Bbb N^{2^r}\leftrightarrow\Bbb N$ that is of degree exactly $2^r$ in all variables, by putting:


*

*$h_0(x)=x$

*$h_{r+1}(x_1,\dots,x_{2^{r+1}})=f(h_r(x_1,\dots,x_{2^r}), h_r(x_{2^r+1},\dots,x_{2^{r+1}}))$

