In Differential Equations and Dynamical Systems by Lawrence Perko, higher-order derivatives are defined in the following way,
Let $E$ be an open subset of $\mathbb{R}^n$. Then for $(\mathbf{x},\mathbf{y})$, we have $D^2\mathbf{f}(\mathbf{x}_0): E \times E \rightarrow \mathbb{R}^n$ is \begin{equation} D^2\mathbf{f}(\mathbf{x}_0)(\mathbf{x},\mathbf{y}) = \sum_{j_1,j_2=1}^n \frac{\partial^2\mathbf{f}(\mathbf{x}_0)}{\partial x_{j_1} \partial x_{j_2}}x_{j_1}y_{j_2}\label{eq:1}\end{equation}
Here is an example given,
Let $\mathbf{f}(\mathbf{x}) = \begin{bmatrix} x_1 + x_1x_2^2 \newline -x_2 + x_2^2 + x_1^2 \end{bmatrix}$. Then, $D^2\mathbf{f}(0,1)(\mathbf{x},\mathbf{y})= \begin{bmatrix} 2x_1y_2 + 2x_2y_1 \newline 2x_1y_1 + 2x_2y_2 \end{bmatrix}$.
I am having trouble understanding how to obtain this result from the given definition. I am presuming the sum in the definition gives the $j_1$ row and $j_2$ column entry. But shouldn't we get a $2 \times 2$ matrix in the example if that was the case?
Thank you in advance!
First of all, your function is vector valued, it has two components $f=(f_1,f_2)$. Your second differential will have two components as well. $D^2f$ is a function of the base point $(x_1,x_2)$ and the increments $(h_1,h_2)$, quadratic in the latter. The notation used in this book is a bit confusing. I would rather write (for a function of two variables) $$ D^2f(x_1,x_2)(h_1,h_2)=\sum_{i,j}\frac{\partial^2f}{\partial x_i\partial x_j}h_ih_j $$ In your case, the Hessian matrix of $f_1$ is $$ H_1=\left( \begin{array}{cc} 0\quad &2x_2\\ 2x_2\quad &2x_1 \end{array} \right) $$ and therefore $$ H_1(0,1)=\left( \begin{array}{cc} 0\quad &2\\ 2\quad &0 \end{array} \right) $$ and then $$ D^2f_1(0,1)(h_1,h_2)=2h_1h_2+2h_2h_1. $$
That is what the first component in your answer means. Same for the second component. BTW, I think there is a typo in the second component of $f$.