How correctly apply a Taylor expansion of first order to a multivariate function decomposition. Suppose I have two function $f(u,v):\mathbb{R}^2\to\mathbb{R}$ and $g(r,s):\mathbb{R}^2\to\mathbb{R}$. I'm interested in the first order Taylor expansion of
$$h:=f(u,g(r,s))$$
To be precise, for a given points $(u_1,v_1)$ and $(r_1,s_1)$ I'm interested in the expansion around $(u_1,g(r_1,s_1))$, i.e. $g(r_1,s_1)= v_1$, evaluated at $(u_1 +\delta (u_2),g(r_1+\delta (r_2),s_1+\delta (s_2)))$, where $\delta( u_2)=u_2-u_1$. The subscript can be seen as a time variable denoting a state. 
Now I want to apply Taylor, but I'm not entirely sure if it is correct. I want to expand the difference to get rid of the constant term:
$$f(u_1+\delta(u_2),g(r_1+\delta (r_2),s_1+\delta (s_2)))-f(u_1,g(r_1,s_1))=\frac{\partial f(u,v)}{\partial u}\bigg|_{=(u_1,g(r_1,s_1))}\cdot[u_2-u_1]+\frac{\partial f(u,v)}{\partial r}\bigg|_{(u_1,g(r_1,s_1))}\cdot\left[g(r_1+\delta (r_2),s_1+\delta (s_2))-g(r_1,s_1)\right]+\frac{\partial f(u,v)}{\partial s}\bigg|_{(u_1,g(r_1,s_1))}\cdot\left[g(r_1+\delta (r_2),s_1+\delta (s_2))-g(r_1,s_1)\right]$$
However we need to apply the chain rule to the two last terms:
$$\frac{\partial f(u,v)}{\partial r}\bigg|_{(u_1,g(r_1,s_1))}=\frac{\partial f(u,v)}{\partial v}\bigg|_{=(u_1,g(r_1,s_1))}\cdot\frac{\partial g(r,s)}{\partial r}\bigg|_{=g(r_1,s_1)}$$
and similar for the other term. About the latter I'm unsure if I evaluate the functions at the right point and if everything is correct. If we use the notation $\delta(g):=g(r_1+\delta (r_2),s_1+\delta (s_2))-g(r_1,s_1)$ we would end with
$$\frac{\partial f(u,v)}{\partial u}\bigg|_{=(u_1,g(r_1,s_1))}\cdot[u_2-u_1]+\left[\frac{\partial g(r,s)}{\partial s}\bigg|_{=g(r_1,s_1)}+\frac{\partial g(r,s)}{\partial r}\bigg|_{=g(r_1,s_1)}\right]\frac{\partial f(u,v)}{\partial v}\bigg|_{=(u_1,g(r_1,s_1))}\cdot\delta(g)$$
It would be very helpful if someone could help me.
 A: I think you're putting the cart before the horse, expanding and evaluating before you find the differentials. We have $h:\mathbb{R}^3\to\mathbb{R}$ defined by $h(x,y,z) = f(x,g(y,z))$ where (following your notation) $f = f(u,v)$ and $g = g(r,s)$. By the multivariable chain rule, the differential will be
$$\frac{\partial h}{\partial x} = \frac{\partial f}{\partial u},$$
$$\frac{\partial h}{\partial y} = \frac{\partial f}{\partial v}\frac{\partial g}{\partial r}$$
$$\frac{\partial h}{\partial z} = \frac{\partial f}{\partial v}\frac{\partial g}{\partial s}.$$
When evaluating at a point, $f$ will be evaluated at $(u,v) = (u,g(r,s))$ and $g$ will be evaluated at $(r,s)$. If the point is $(x_0,y_0,z_0)$, the second two coordinates are plugged into $g$, so the $f$ derivatives will be evaluated at $(x_0,g(y_0,z_0))$ and the $g$ derivatives will be evaluated at $(y_0,z_0)$.
So the first order Taylor expansion at $(x_0,y_0,z_0)$ will be
\begin{align*}
 f(x_0,g(y_0,z_0)) &+ \frac{\partial f}{\partial u}(x_0,g(y_0,z_0))(x-x_0)\\
 &+ \frac{\partial f}{\partial v}(x_0,g(y_0,z_0))\frac{\partial g}{\partial r}(y_0,z_0)(y-y_0)\\
&+ \frac{\partial f}{\partial v}(x_0,g(y_0,z_0))\frac{\partial g}{\partial s}(y_0,z_0)(z-z_0)
\end{align*}
