Understanding second derivatives 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
 A: 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}$$
