Deriving 2-variable Taylor expression from 1-variable Suppose we know Taylor's theorem for $f: \mathbb{R} \rightarrow \mathbb{R}$ with remainder. How would we use this to derive Taylor's theorem for $g: \mathbb{R}^2 \rightarrow \mathbb{R}$ with remainder?
To start, I'm trying to figure out which version of the 1-dimensional Taylor remainder to use (options are Lagrange, Cauchy, or integral form). To do this, I need to figure out the remainder for 2-variable Taylor. Wikipedia lists the remainder for a general n-variable Taylor's theorem, but I'm having trouble extracting the 2-dimensional case from it.
 A: (this is pretty much what Michael said, I just happen to have a few notes on this and I thought this might be helpful. If I were you I would try to extrapolate from Lagrange's form of the remainder)
We already know Taylor's theorem for functions on $\mathbb{R}$,
$$ g(x) = g(a)+g'(a)(x-a) + \frac{1}{2}g''(a)(x-a)^2 +\cdots +\frac{1}{k!}g^{(k)}(a)(x-a)^k + R_k $$
and... If the remainder term vanishes as $k \rightarrow \infty$ then the function $g$ is represented by the Taylor series given above and we write:
$$ g(x) = \sum_{k=0}^{\infty}\frac{1}{k!}g^{(k)}(a)(x-a)^k. $$
Consider the function of two variables $f: U \subseteq \mathbb{R}^2 \rightarrow \mathbb{R}$ which is smooth with smooth partial derivatives of all orders. Furthermore, let $(a,b) \in U$ and construct a line through $(a,b)$ with direction vector $(h_1,h_2)$ as usual:
$$ \phi (t) = (a,b) + t(h_1,h_2) = (a+th_1, b+th_2) $$
for $t \in \mathbb{R}$. Note $\phi(0)=(a,b)$ and $\phi '(t) = (h_1,h_2) = \phi ' (0)$. Construct $g = f \,{\scriptstyle \stackrel{\circ}{}}\, \phi  : \mathbb{R} \rightarrow \mathbb{R}$ and choose $dom(g)$ such that $\phi(t) \in U$ for $t \in dom(g)$. This function $g$ is a real-valued function of a real variable and we will be able to apply Taylor's theorem from calculus II on $g$. However, to differentiate $g$ we'll need tools from calculus III to sort out the derivatives. In particular, as we differentiate $g$, note we use the chain rule for functions of several variables:
$$
\begin{align} \notag
   g'(t) =(f \,{\scriptstyle \stackrel{\circ}{}}\, \phi)'(t) &= f'(\phi(t))\phi '(t)  \\ \notag
         &= \nabla f(\phi(t)) \cdot (h_1,h_2) \\ \notag
         &= h_1f_x(a+th_1, b+th_2)+h_2f_y(a+th_1, b+th_2) \notag
\end{align}
$$
Note $g'(0)=h_1f_x(a, b)+h_2f_y(a, b)$. Differentiate again (I omit $(\phi(t))$ dependence in the last steps),
$$
\begin{align} \notag
   g''(t) &= h_1f_x '(a+th_1, b+th_2)+h_2f_y '(a+th_1, b+th_2)  \\ \notag
          &= h_1 \nabla f_x (\phi(t)) \cdot (h_1,h_2) +h_2 \nabla f_y (\phi(t)) \cdot (h_1,h_2)  \\ \notag
          &= h_1^2f_{xx}+ h_1h_2f_{yx} + h_2h_1f_{xy}+ h_2^2f_{yy}  \\ \notag
          &= h_1^2f_{xx}+ 2h_1h_2f_{xy} + h_2^2f_{yy}  \notag
\end{align}
$$
Thus, making explicit the point dependence, $g''(0) = h_1^2f_{xx}(a,b)+ 2h_1h_2f_{xy}(a,b) + h_2^2f_{yy}(a,b)$. We may construct the Taylor series for $g$ up to quadratic terms:
$$
\begin{align} \notag
   g(0+t) &= g(0)+tg'(0)+\frac{1}{2}g''(0) + \cdots \\ \notag
&= f(a,b)+t[h_1f_x(a, b)+h_2f_y(a, b)]+\frac{t^2}{2}\bigl[ h_1^2f_{xx}(a,b)+ 2h_1h_2f_{xy}(a,b) + h_2^2f_{yy}(a,b) \bigr] + \cdots \\ \notag
\end{align}
$$
Note that $g(t)= f(a+th_1, b+th_2)$ hence $g(1) = f(a+h_1, b+h_2)$ and consequently,
$$
\begin{align} \notag
 f(a+h_1, b+h_2) &= f(a,b)+h_1f_x(a, b)+h_2f_y(a, b)+ \\ \notag 
                 & \qquad +\frac{1}{2}\biggl[ h_1^2f_{xx}(a,b)+ 2h_1h_2f_{xy}(a,b) + h_2^2f_{yy}(a,b) \biggr] + \cdots  \notag
\end{align}
$$
Omitting point dependence on the $2^{nd}$ derivatives,
$$ \boxed{ f(a+h_1, b+h_2) = f(a,b)+h_1f_x(a,b)+h_2f_y(a,b)+ 
                 \tfrac{1}{2}\bigl[ h_1^2f_{xx}+ 2h_1h_2f_{xy} + h_2^2f_{yy} \bigr] + \cdots }$$
Sometimes we'd rather have an expansion about $(x,y)$. To obtain that formula simply substitute $x-a = h_1$ and $y-b = h_2$. Note that the point $(a,b)$ is fixed in this discussion so the derivatives are not modified in this substitution,
$$
\begin{align} \notag
 f(x, y) &= f(a,b)+(x-a)f_x(a, b)+(y-b)f_y(a, b)+ \\ \notag 
                 & \qquad +\frac{1}{2}\biggl[ (x-a)^2f_{xx}(a,b)+ 2(x-a)(y-b)f_{xy}(a,b) + (y-b)^2f_{yy}(a,b) \biggr] + \cdots  \notag
\end{align}
$$
At this point we ought to recognize the first three terms give the tangent plane to $z = f(z,y)$ at $(a,b,f(a,b))$. The higher order terms are nonlinear corrections to the linearization, these quadratic terms form a quadratic form. If we computed third, fourth or higher order terms we will find that, using $a=a_1$ and $b=a_2$ as well as $x=x_1$ and $y=x_2$,
$$ \boxed{f(x, y) = \sum_{n=0}^{\infty} \sum_{i_1=0}^{2}\sum_{i_2=0}^{2} \cdots \sum_{i_n=0}^{2} \frac{1}{n!}
\frac{\partial^{(n)}f(a_1,a_2)}{\partial x_{i_1}\partial x_{i_2} \cdots \partial x_{i_n}}
(x_{i_1} -a_{i_1})(x_{i_2} -a_{i_2})\cdots (x_{i_n} -a_{i_n})} $$
A: This probably isn't totally what you want, but I'll leave it here anyway as it may help you on your way.
This proves the second-order Taylor formula for scalar fields. I think the general case just involves some pencil pushing to figure out the $n^\text{th}$ derivative of $f$.
Let $g:\Bbb R^2 \to \Bbb R$. Then consider $f(t)=g(\mathbf a + t \mathbf y)$ for $|t|\le 1$. Then apply the Taylor formula to $f$, i.e. we have that $$f(1)=f(0)+f'(0)+\frac 1 {2!}f''(c), \quad\text{where}\quad 0<c<1$$
(Here we used Lagrange's form for the remainder.) Now since $f$ is a composite function, we have that $f'(t)=\nabla g(\mathbf a + t\mathbf y)\cdot \mathbf y$. In particular, $f'(0)=\nabla g(\mathbf a)\cdot \mathbf y$. Now all that's left is to differentiate $f'(t)$ again. Write
$$\begin{array}{rl}f'(t)=&D_1 g(\mathbf a+ t \mathbf y)y_1+D_2g(\mathbf a +t\mathbf y)y_2\\
\implies f''(t)=&D_1D_1g(\mathbf a + t \mathbf y)y_1^2+D_2D_1g(\mathbf a +t\mathbf y)y_2y_1+\\
&D_1D_2g(\mathbf a + t\mathbf y)y_1y_2+D_2D_2g(\mathbf a + t \mathbf y)y_2^2\\
=&\sum_{i=1}^2\sum_{j=1}^2D_{ij}g(\mathbf a + t \mathbf y)y_iy_j\end{array}$$
Note that we may write $f''(t)=\mathbf y H[\mathbf a + t\mathbf y] \mathbf y^\text T$, where $H$ is the Hessian Matrix of $f$. Thus we have (from the Taylor formula for $f$):
$$g(\mathbf a+\mathbf y)=g(\mathbf a)+\nabla g(\mathbf a )\cdot \mathbf y + \frac 1 2 \mathbf y H[\mathbf a+c\mathbf y]\mathbf y^\text T$$
It seems like this can be extended further, by a general formula for the calculation of the $n^\text{th}$ derivative of $f$, which is what they do on the Wikipedia page. Also note that here I used Lagrange's form of the remainder, while Wikipedia has the Integral form of the remainder. Essentially any form of the remainder will work, however certain forms may be more useful than others for particular applications.
