I am having a hard time understanding how to determine the second derivative of a matrix.
I have researched Hessian matrices and do not see how i would apply it to vector funciton.
problem statement: compute the derivative of the following
$$ f(x) =\begin{bmatrix} x_1+x_1x_2^2\\ -x_2+x_2^2+x_1^2\\ \end{bmatrix}$$
I have found
$$ DF =\begin{bmatrix} 1+x_2^2&2x_1\\ 2x_1&-1+x_2\\ \end{bmatrix} $$
I now need to compute $D^2f(x_0)(x,y)$. I have
$$ D^2f(x_0)(x,y) = \sum_{j_1,j_2=1}^n \frac{\partial^2f(x_0)}{\partial x_{j_1} \partial x_{j_2}} x_{j_1}y_{j_2} $$
But don't quite understand it. I have also searched Hessian Matrix, however the majority of times when it is used is when there is only one equation.
Could someone please walk me through the $$D^2f(x_0)(x,y)$$ equation and also explain if a Hessian matrix is applicable here.
I will provide $x_0$ if needed
Note: I'm assuming your question is about how to calculate the second derivative as opposed to what it is. Thus I will not go cover the isomorphism $\mathcal L(X,\mathcal L(X,Y)) \simeq \operatorname{Bil}(X\times X\to Y)$ nor where the formula for calculating the derivative in Cartesian coordinates comes from.
The second derivative of a vector field $f:\Bbb R^m \to \Bbb R^n$ at a point $\mathbf a$ is a bilinear function $D^2f(\mathbf a):\Bbb R^m\times \Bbb R^m \to \Bbb R^n$. If it exists, then in Cartesian coordinates it is given by $$D^2f(\mathbf a)(\mathbf u,\mathbf v) = \sum_{j_1,j_2=1}^n \frac{\partial^2f(\mathbf a)}{\partial x_{j_2} \partial x_{j_1}} u_{j_1}v_{j_2}$$ as you wrote.
We'll use your function as an example. Let $f:\Bbb R^2 \to \Bbb R^2$ be given by $f(\mathbf x) = (x+xy^2,-y+y^2+x^2)$, where $\mathbf x = (x,y)$. Then $$\partial_x f(\mathbf x) = (1+y^2,2x) \\ \partial_yf(\mathbf x) = (2xy,-1+2y)$$ $$\begin{matrix}\partial_{xx} f(\mathbf x)= (0,2) & \partial_{yx} f(\mathbf x) = (2y,0) \\ \partial_{xy} f(\mathbf x) = (2y,0) & \partial_{yy}f(\mathbf x) = (2x,2)\end{matrix}$$
So then the second derivative at $\mathbf a = (6,1)$ evaluated at, say, $(\mathbf u,\mathbf v)$ where $\mathbf u = (2,1)$ and $\mathbf v = (\frac 13,\frac 12)$ is
$$\begin{align}D^2f(6,1)(\mathbf u,\mathbf v) &= \partial_{xx}f(6,1)u_xv_x + \partial_{xy}f(6,1)u_yv_x + \partial_{yx}f(6,1)u_xv_y + \partial_{yy}f(6,1)u_yv_y \\ &= (0,2)\cdot 2 \cdot \frac 13 + (2,0)\cdot 1\cdot \frac 13 + (2,0)\cdot 2 \cdot \frac 12 + (12,2)\cdot 1\cdot \frac 12 \\ &= \left(0,\frac 43\right) + \left(\frac 23,0\right) + (2,0) + (6,1) \\ &= \frac 13(26,7)\end{align}$$