Let $\mathbf{\vec{x}} = \left[\begin{array}{c} x_1 \\ x_2 \\ x_3 \\ . \\ . \\ . \\ x_n \end{array}\right]$ and $f(\mathbf{\vec{x}}): \mathbb{R}^n \to \mathbb{R}$, $n \in \mathbb{N}$ be a scalar-valued function, then its Quadratic approximation $Q(\mathbf{\vec{x}})$ about a point $\mathbf{\vec{x}}_0$ is given by:
$$Q(\mathbf{\vec{x}}) = f(\mathbf{\vec{x}}_0) + \mathbf{\nabla}f(\mathbf{\vec{x}}_0)^T (\mathbf{\vec{x} - \vec{x}}_0) + \dfrac{1}{2} (\mathbf{\vec{x} - \vec{x}}_0)^T \mathbf{H}_f (\mathbf{\vec{x}}_0) (\mathbf{\vec{x} - \vec{x}}_0)$$
where $\mathbf{\nabla}f$ is the gradient and $\mathbf{H}_f$ the Hessian matrix of $f$.
We know that the function $Q(\mathbf{\vec{x}})$ is suppose to have the same output, first and second partial derivatives as $f(\mathbf{\vec{x}})$ at the point $\mathbf{\vec{x}}_0$.
Clearly, $Q(\mathbf{\vec{x}}_0) = f(\mathbf{\vec{x}}_0)$. Now, how do I prove the following?:
$$\mathbf{\nabla}Q(\mathbf{\vec{x}}_0) = \mathbf{\nabla}f(\mathbf{\vec{x}}_0)$$
$$\mathbf{H}_Q(\mathbf{\vec{x}}_0) = \mathbf{H}_f(\mathbf{\vec{x}}_0)$$
You can verify that the gradients and Hessian agree by direct calculation. Observe that both $\nabla f[\vec{\mathbf x}_0]$ and ${\mathbf H}_f[\vec{\mathbf x}_0]$ are constants, so the second term of $Q$ is a linear function of $\vec{\mathbf x}-\vec{\mathbf x}_0$ and the third term a quadratic form in $\vec{\mathbf x}-\vec{\mathbf x}_0$. (I use square brackets to indicate evaluation at a point in an attempt at disambiguation in the following.)
Analogous to the single-variable differentiation rule $(cx)'=c$ we have $$\nabla(\vec{\mathbf v}^T\vec{\mathbf x}) = \nabla\left(\sum_i v_ix_i\right) = \vec{\mathbf v},$$ therefore $$\nabla(\nabla f[\vec{\mathbf x}_0]^T(\vec{\mathbf x}-\vec{\mathbf x}_0)) = \nabla(\nabla f[\vec{\mathbf x}_0]^T\vec{\mathbf x})-\nabla(\nabla f[\vec{\mathbf x}_0]^T\vec{\mathbf x}_0) = \nabla f[\vec{\mathbf x}_0].$$ More generally, the differential of a linear function $A\mathbf x$ of $\mathbf x$ is $A$—the best linear approximation to a linear map is the map itself.
For the quadratic term, we can use the product rule $\nabla(u^Tv)=u^T\nabla v+(\nabla u)^Tv$: $$\begin{align} \nabla\left(\frac12(\vec{\mathbf x}-\vec{\mathbf x}_0)^T{\mathbf H}_f[\vec{\mathbf x}_0](\vec{\mathbf x}-\vec{\mathbf x}_0)\right) &= \frac12 (\vec{\mathbf x}-\vec{\mathbf x}_0)^T \nabla({\mathbf H}_f[\vec{\mathbf x}_0](\vec{\mathbf x}-\vec{\mathbf x}_0)) + \frac12(\nabla(\vec{\mathbf x}-\vec{\mathbf x}_0))^T{\mathbf H}_f[\vec{\mathbf x}_0](\vec{\mathbf x}-\vec{\mathbf x}_0) \\ &= \frac12 (\vec{\mathbf x}-\vec{\mathbf x}_0)^T {\mathbf H}_f[\vec{\mathbf x}_0] + \frac12 {\mathbf H}_f[\vec{\mathbf x}_0](\vec{\mathbf x}-\vec{\mathbf x}_0) \\ &= {\mathbf H}_f[\vec{\mathbf x}_0](\vec{\mathbf x}-\vec{\mathbf x}_0) \end{align}$$ since ${\mathbf H}_f$ is symmetric. This vanishes at $\vec{\mathbf x}_0$, therefore $\nabla Q[\vec{\mathbf x}_0] = \nabla f[\vec{\mathbf x}_0]$. (We could’ve gotten to this without actually differentiating by noting that the result must be a linear combination of $x_i-x_{0i}$ terms, all of which vanish.)
For ${\mathbf H}_Q[\vec{\mathbf x}_0]$ we can proceed in a similar fashion by differentiating $\nabla Q = \nabla f[\vec{\mathbf x}_0]+{\mathbf H}_f[\vec{\mathbf x}_0](\vec{\mathbf x}-\vec{\mathbf x}_0)$ directly, producing ${\mathbf H}_f[\vec{\mathbf x}_0]$ per the earlier observation regarding the differential of a linear function, or note that the Hessian of a constant or linear function vanishes, while the Hessian of a quadratic form is 2 times the (constant) matrix of the form (analogously to the single-variable $(ax^2)''=2a$), so that the only term that survives when computing the Hessian directly from $Q$ is ${\mathbf H}_f[\vec{\mathbf x}_0]$.