How to easily create a polynomial function that gives a desired output? I am looking for an easy way (formula or algorithm) to create a polynomial function that gives the arbitrary preset output for the first values of x. For instance, the desired output can be $y = 1, 2, 3, 4, 5, 6, 100$ for $x = 1, 2, 3, 4, 5, 6, 7$. I can create one of the many functions by hand as: 

$y = x + (x-1)(x-2)(x-3)(x-4)(x-5)(x-6)\cdot \frac{93}{720}.$

But for other desired outputs it is really painstaking. For example, I am working now with the desired output $y = 2, 10, 12, 16, 17, 18, 19, 200, 300$ when $x = 1, 2, 3, 4, 5, 6, 7, 8, 9$.
Notes:  
$x$ is always a natural positive number ($x > 0$). It doesn't matter what happens for other values of $x$ (greater than the ones given).
- I am looking for a easy way, as I don't know / I don't have / I have no experience with math related software.
- Also for a polynomial (or similar simple) function, as my mathematical background is quite limited, so I don't know any gamma, delta or other functions, series, etc.
- This quest is simply for recreation, to show friends how there are different solutions to the problems of the type: "Find the next term in this succession" 
 A: Polynomials in the standard form (as a sum of monomials) are actually the 'wrong' sort of 'encoding' to do interpolation. This post explains how to very efficiently interpolate any polynomial function, without ever computing the polynomial coefficients! In short, we can easily express any given polynomial sequence as a Newton series by computing successive forward differences of the given sequence, and not only can we use this to efficiently compute subsequent terms in the sequence, we can also quickly convert this representation back to standard polynomial form, in the same way the linked post demonstrates for the sum of the first $n$ perfect cubes.
For your example the iterated forward differences are easily computed by hand:
 1,  2,  3,  4,  5,  6,100,659,...
 1,  1,  1,  1,  1, 94,559,...
 0,  0,  0,  0, 93,465,...
 0,  0,  0, 93,372,...
 0,  0, 93,279,...
 0, 93,186,...
93, 93,...

Which also means that you can use this forward difference technique (which works for any polynomial sequence such as $1,4,9,16,25,\cdots$) to justify that you can extend any given sequence polynomially as a reasonable pattern.
That said, if too few terms of a sequence are given (which is often the case), it is actually easier to describe it as simply a list with no pattern, because the intended pattern would have a description that is longer than just listing the few terms given.
A: As mentioned in a comment, the key is "Lagrange Interpolation". (See, for instance, Wikipedia or MathWorld.) The formal description can be a little opaque, but the idea is actually pretty straightforward: Build a polynomial out of parts, most of which go away with each input, then make the remaining part do what it needs to do.
For example, suppose you want the three (distinct) inputs $a$, $b$, $c$ to give the three outputs $p$, $q$, $r$. Lagrange claims that the polynomial you want has the form
$$h(x-b)(x-c) \;+\; j(x-c)(x-a) \;+\; k(x-a)(x-b) \tag{$\star$}$$
where $h$, $j$, $k$ are some as-yet-undetermined constants. 
When $x=a$, two terms of $(\star)$ vanish, leaving $h(a-b)(a-c)$; likewise, $x=b$ yields $j(b-c)(b-a)$, while $x=c$ yields $k(c-a)(c-b)$. So, just choose $h$, $j$, $k$ to give the outputs you desire:
$$\begin{align}
h(a-b)(a-c) = p \qquad\to\qquad h = \frac{p}{(a-b)(a-c)} \\[4pt]
j(b-c)(b-a) = q \qquad\to\qquad \,j = \frac{q}{(b-c)(b-a)} \\[4pt]
k(c-a)(c-b) = r \qquad\to\qquad k = \frac{r}{(c-a)(c-b)}
\end{align}$$
Thus, we can re-write $(\star)$ informally as
$\require{cancel}$
$$p\;\frac{\cancel{(x-a)}(x-b)(x-c)}{\cancel{(a-a)}(a-b)(a-c)}
\;+\;q\;\frac{(x-a)\cancel{(x-b)}(x-c)}{(b-a)\cancel{(b-b)}(b-c)}
\;+\;r\;\frac{(x-a)(x-b)\cancel{(x-c)}}{(c-a)(c-b)\cancel{(c-c)}} \tag{$\star\star$}$$
Clearly, this kind of thing can be extended to any number of inputs and outputs.
A: As Gabriel Romon commented: Lagrange interpolation does the trick. 
The polynomials form a vector space (which works very similar to the 3D coordinate space) where the coefficients work as coordinates and the monomials including $x^0=1$ as orthogonal standard base. The Lagrange polynomials form another one. They are constructed from your given sampling points ($x$-values). The nice thing is that each Lagrange polynomial is zero in all but one sampling points. And in this sampling point it has the value $1$. So you can multiply each of the Lagrange polynomials with the function value at the sampling point where it has a nonzero value and sum over all these polynomials.
This looks as easy as 
$$
y = \sum_{i=1}^n y_i \cdot \ell_i(x).
$$
As you can see there is no $x_i$ in this formula any more. This makes it possible to precompute the Lagrange polynomials $\ell_i$ and solve all your problems together in a secound step, as long as all functions are defined using the same sampling points.
