Yes, you can certainly just do the second order Taylor expansions of $f_1$ and $f_2$ and then just wrap each like term together in a matrix:
$$f_1(x+h_1,y+h_2) \approx f_1(x,y) + [Df_1(x,y)](h_1,h_2) + \frac 12[D^2f_1(x,y)](h_1,h_2,h_1,h_2) \\ f_2(x+h_1,y+h_2) \approx f_2(x,y) + [Df_2(x,y)](h_1,h_2) + \frac 12[D^2f_2(x,y)](h_1,h_2,h_1,h_2) \\ \implies \pmatrix{f_1(x+h_1, y+h_2) \\ f_2(x+h_1, y+h_2)} \approx \pmatrix{f_1(x,y) \\ f_2(x,y)} + \pmatrix{[Df_1(x,y)](h_1,h_2) \\ [Df_2(x,y)](h_1,h_2)} + \frac 12\pmatrix{[D^2f_1(x,y)](h_1,h_2,h_1,h_2) \\ [D^2f_2(x,y)](h_1,h_2,h_1,h_2) }$$
where, in my notation $$[Df_i(x,y)](h_1,h_2) = \nabla f_i(x,y)\cdot (h_1,h_2) \\ [D^2f_i(x,y)](h_1,h_2,h_1,h_2) = \pmatrix{h_1 &h_2}[Hf_i(x,y)]\pmatrix{h_1 \\ h_2}$$ and $Hf_i(x,y)$ is the Hessian matrix of $f_i$ at $(x,y)$.
But notice that this is just
$$f(x+h_1, y+h_2) \approx f(x,y) + [Jf(x,y)](h_1,h_2) + \frac 12\color{red}{\pmatrix{[D^2f_1(x,y)](h_1,h_2,h_1,h_2) \\ [D^2f_2(x,y)](h_1,h_2,h_1,h_2) }}$$ where $Jf$ is the Jacobian of $f$ at $(x,y)$ and where you've probably never seen that last term before.
I chose to write the second derivative (evaluated at $(h_1,h_2,h_1,h_2)$) in terms of its components because I suspect you haven't gotten to higher than first derivatives of functions between normed spaces. Higher derivatives of a function $f: \Bbb R^m \to \Bbb R^n$ are really mappings of the type $D^2f: \mathbb{R}^m \rightarrow \mathcal{L}(\mathbb{R}^m,\mathcal{L}(\mathbb{R}^m,\mathbb{R}^n))$. However, after we choose a point to evaluate the second derivative at we can use currying to associate this with a bilinear function of type $\Bbb R^m \times \Bbb R^m \to \Bbb R^n$. Unfortunately there is no way to multiply a matrix by two column matrices to get another column matrix with the correct components. What we need here, if we want to condense the notation for that last object, are tensors.
Think of the Taylor series expansion of $f(x) = e^x$ around $x=a$.
$f(x) = f(a) + f'(a)(x-a) + {f''(a) \over 2!}(x-a)^2 + \ldots $
Using the fact that all derivatives of $f(x) = e^x$ are also $e^x$, we get:
$e^x = f(x) = e^a + e^a(x-a) + {e^a \over 2!}(x-a)^2 + \ldots$
Setting $a=0$, $e^x = 1 + x + {x^2 \over 2!} + \ldots$ (Maclaurin series)
Now, let's consider the MGF you computed:
$m(t) = e^{-\lambda} \sum_y \frac{(\lambda e^t)^y}{y!} $
Set $x = \lambda e^t$ to get:
$m(t) = e^{-\lambda} \sum_y {x^y \over y!} = e^{-\lambda} \left( 1 + x + {x^2 \over 2!} + \ldots \right) = e^{-\lambda} e^x = e^{-\lambda} e^{\lambda e^t}$ (using our Maclaurin series)
Best Answer
Yes, the product of those series equals $e^{-(x^2+y^2)}.$ But that product is not quite the Taylor series of this function centered at $(0,0);$ that would be
$$\sum_{m=0}^{\infty}\sum_{n=0}^{\infty} (-1)^{m+n}\frac{x^{2m}}{m!}\frac{y^{2n}}{n!}.$$
The same ideas apply at $(1,2).$