How to do a Taylor expansion of a vector-valued function 
Let $f:\Bbb R^2\to \Bbb R^2$ be given by $$f(x,y):= \left(e^x\sin(x+y),e^{y-x}\tanh(y)\right)$$ Find the second-order Taylor expansion of $f$ about (x,y)=(0,0)$.

I know how to find the Taylor Expansion of a function $f:\Bbb R^n \to \Bbb R$. But I have never seen $f:\Bbb R^n \to \Bbb R^m$ before.  Do I find the Taylor expansion of
$$f_1(x,y)= \sin(x+y)e^{x}\quad  \text{at } (0,0)$$
and 
$$f_2(x,y) = \tanh(y)e^{y-x}\quad \text{at } (0,0)$$
and then combine the expansions for the Taylor expansion of $f$?
 A: 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.
