Taylor expansion involving vectors 
Question: 
Consider a series of functions $f_1,\cdots,f_d:\mathbb{R}^d\to\mathbb{R}$ such that \begin{align}f_1(v_1,\cdots,v_d)&=0\\
&\vdots\\
f_d(v_1,\cdots,v_d)&=0\end{align}
for $v_1,v_2,\cdots,v_d\in \mathbb{R}$
Fix $\textbf{x},\textbf{y}\in \mathbb{R}^d$ and let $$g(\lambda)=f_i(\textbf{x}+\lambda(\textbf{y}-\textbf{x}))$$ for $\lambda \in \mathbb{R}$
Write down a Taylor series for $g$ near $\lambda=0$ to second order.

I know that this Taylor expansion will have the form $$g(\lambda)=g(0)+\lambda g'(0)+\frac{\lambda ^2}{2}g''(0)$$ (differentiating wrt $\lambda$) however I am struggling to actually compute this.
I have computed the following:
\begin{align}g'(\lambda)&=(\textbf{y}-\textbf{x}) f_i '(\textbf{x}+\lambda (\textbf{y}-\textbf{x}))\\
g''(\lambda)&=(\textbf{y}-\textbf{x})^2 f_i ''(\textbf{x}+\lambda (\textbf{y}-\textbf{x}))\end{align}
however $g''(\lambda)$ involves the square of a vector, which I feel is incorrect.
Can someone either point me in the right direction or explain to me where I am going wrong please.

Edit: I found this related question and have come up with a potential solution. Can someone verify it for me please?
$$g(\lambda) = f_i(\textbf{x}) + \lambda \sum_{j=1}^d \frac{\partial f_i(\textbf{x})}{\partial x_j}(y_j-x_j) + \frac{\lambda^2}2 \sum_{j,k=1}^d\frac{\partial f_i(\textbf{x})}{\partial x_j\partial x_k}(y_j-x_j)(y_k-x_k)$$
 A: These expansions get confusing quickly with vector valued variables and functions as arguments, so I'll renotate a bit for (hopefully) more clarity and derive the result using just the multivariate chain rule.
Let $\lambda\in\mathbb{R}$, $f_i : \mathbb{R}^n \rightarrow\mathbb{R}$, $u : \mathbb{R} \rightarrow\mathbb{R}^n$, $x,y\in\mathbb{R}^n$, and 
$$ g(\lambda)=f_i(x+\lambda (y-x)) = f_i(u(\lambda))$$
where $u(\lambda)=x+\lambda(y-x)$.
Notice that $\partial_\lambda u = y-x$.
Then the second-order univariate Taylor expansion about $\lambda=a$ is $$ T_g(\lambda|a) = g(a)+[\lambda - a]\partial_\lambda g(a) + \frac{1}{2}\partial_{\lambda\lambda} g(a) [\lambda - a]^2 $$
where $\partial_\lambda=\partial/\partial \lambda$. We really just need $\partial_{\lambda} g(a)$ and $\partial_{\lambda\lambda} g(a)$.
The first one, using the multivariate chain rule:
\begin{align}
\partial_\lambda g(\lambda) 
&= \sum_k \left.\frac{\partial f_i}{\partial u_k} \right\rvert_{u(\lambda)} \partial_\lambda u_k(\lambda) \\
&= \nabla f_i(u(\lambda))^T [y-x] \tag{1}
\end{align}
Now from (1) we can evaluate at the special case of $\lambda=0$:
\begin{align}
\partial_\lambda g(0) &= \sum_k \left.\frac{\partial f_i}{\partial u_k} \right\rvert_{u(0)} \partial_\lambda u_k(0) = \sum_k \left. \frac{\partial f_i}{\partial x_k}\right\rvert_{x} [y_k - x_k] \\ &= \sum_k  [y_k - x_k] \frac{\partial}{\partial x_k}  f_i(x) = (y-x)^T \nabla_x f_i(x) \tag{2} 
\end{align}
because $$ \left.\frac{\partial u_k}{\partial x_k}\right|_{\lambda=0} = \left.(1-\lambda)\right|_{\lambda=0} =1 $$ so
$$ \left.\frac{\partial f_i}{\partial x_k}\right|_{\lambda=0} = \left.\frac{\partial f_i}{\partial u_k}\right|_{u(0)} \left. \frac{\partial u_k}{\partial x_k}\right|_{\lambda = 0} \implies \left.\frac{\partial f_i}{\partial u_k}\right|_{u(0)} = \left.\frac{\partial f_i}{\partial x_k}\right|_{x}  $$
Ok, now for the second derivative. First, I want to define $q :\mathbb{R}^n \rightarrow \mathbb{R}^n$ as $q(u(\lambda)) = \nabla f_i(u(\lambda))$, meaning that $q_j : \mathbb{R}^n \rightarrow \mathbb{R}$ can be written $q_j(u(\lambda)) = (\partial f_i/\partial u_j)|_{u(\lambda)}$. This means that the Hessian matrix $\mathcal{H}$ of $f$ has rows given by $$\nabla q_j(u(\lambda))^T = \left. \left[ \frac{\partial^2 f_i}{\partial u_1 \partial u_j},\cdots,\frac{\partial^2 f_i}{\partial u_n \partial u_j} \right]\right\rvert_{u(\lambda)}, $$
that is, $ \mathcal{H}[f_i(u(\lambda))] $ has components given by: 
$$ \mathcal{H}_{\ell k} = \left.\frac{\partial^2 f_i}{\partial u_\ell \partial u_k}\right|_{u(\lambda)} = \nabla q_\ell(u(\lambda))_k $$
Anyway, we can then compute (using the chain rule again):
$$ \frac{\partial}{\partial \lambda} q_j(u(\lambda)) = \nabla q_j(u(\lambda))^T \frac{\partial u}{\partial \lambda} = \nabla q_j(u(\lambda))^T[y-x] $$
Alright, so for the second derivative, start with the result from (1):
\begin{align}
\partial_{\lambda\lambda} g(\lambda)
&= \partial_\lambda \left[ \nabla f_i(u(\lambda))^T [y-x] \right] \\
&= \sum_j (y_j - x_j)\frac{\partial}{\partial\lambda}\frac{\partial}{\partial u_j}f_i(u(\lambda)) \\
&= \sum_j (y_j - x_j)\frac{\partial}{\partial\lambda}q_j(u(\lambda)) \\ 
&= \sum_j (y_j - x_j)\nabla q_j(u(\lambda))^T[y-x] \\ 
&= \sum_j (y_j - x_j) \sum_k \nabla q_j(u(\lambda))_k[y_k-x_k] \\ 
&= \sum_j (y_j - x_j) \sum_k 
\frac{\partial^2  f_i(u(\lambda))}{\partial u_j \partial u_k} [y_k-x_k] \\ 
&= (y-x)^T \mathcal{H}[f_i(u(\lambda))](y-x) \tag{3}
\end{align}
Next, let's use the following relation, which will help us get rid of $u$:
\begin{align}
\frac{\partial^2 f_i}{\partial x_\ell \partial x_k} 
&= \frac{\partial}{\partial x_\ell} \frac{\partial f_i}{\partial x_k}
=  \frac{\partial}{\partial x_\ell} 
   \left(\frac{\partial f_i}{\partial u_k}\frac{\partial u_k}{\partial x_k}\right)
= \underbrace{\frac{\partial u_k}{\partial x_k}}_{1-\lambda} \;
 \underbrace{ \frac{\partial^2 f_i}{\partial x_\ell \partial u_k} }_{\partial q_k(u(\lambda))/ \partial {x_\ell}}
  + 
  \frac{\partial f_i}{\partial u_k}\underbrace{\frac{\partial}{\partial x_\ell}\frac{\partial u_k}{\partial x_k}}_0 \\
&= (1-\lambda)\frac{\partial}{\partial u_k} \frac{\partial f_i}{\partial{x_\ell}}
= (1-\lambda)\frac{\partial}{\partial u_k}\left[
\frac{\partial f_i}{\partial{u_\ell}}\frac{\partial u_\ell}{\partial{x_\ell}}
\right] = (1-\lambda)^2\frac{\partial^2  f_i(u(\lambda))}{\partial u_j \partial u_k}
\end{align}
So, at $\lambda = 0$ we have
$$ \left.\frac{\partial^2  f_i(u(\lambda))}{\partial u_j \partial u_k}\right|_{u(0)} 
= \frac{\partial^2 f_i(x)}{\partial x_\ell \partial x_k} \tag{4} $$
which means that (using (3) and (4))
$$ \partial_{\lambda\lambda} g(0) = \sum_{j,k} (y_j - x_j)  
\frac{\partial^2  f_i(u(0))}{\partial u_j \partial u_k} [y_k-x_k]
= \sum_{j,k} [y_j - x_j] \frac{\partial^2 f_i(x)}{\partial x_\ell \partial x_k} [y_k-x_k] \tag{5} $$
Altogether, using (1) and (3), the general form of the second order Taylor expansion of $g(\lambda)$ expanded about $a$ can be written as
$$
T_g(\lambda|a) = g(a) + [\lambda - a]\nabla f_i(u(a))^T [y-x] + \frac{1}{2}[\lambda - a]^2 (y-x)^T \mathcal{H}[f_i(u(a))](y-x)
$$
and in the desired special case where $a=0$ we get (using (2) and (5)):
\begin{align}
T_g(\lambda|0) 
&= f_i(x) + \lambda (y-x)^T \nabla_x f_i(x) + \frac{\lambda^2}{2} (y-x)^T \mathcal{H}[f_i(x)](y-x) \\
&= f_i(x) + \lambda \sum_\ell (y_\ell - x_\ell)\partial_x f_i(x) + \frac{\lambda^2}{2}
\sum_{j,k}(y_j-x_j)(y_k - x_k) \frac{\partial^2 f_i(x)}{\partial x_j \partial x_k}
\end{align}
exactly as you did. :D
