I'm trying to implement the algorithm below in matlab but i dont know how to compute the gradient and hessian of a vector valued function, i need to know how to do it to apply what is in the red underlined parts.

Lets say $h(x) = ( (x_{1} +1)^2 + x_{2}^2 - 1, x_{1}^2 + x_{2}^2 -2, x_{1}^3) $ , how can I find the expression for the gradient and hessian of $h$?
$ \def\n{\nabla} \def\LR#1{\left(#1\right)} \def\m#1{\left[\begin{array}{r}#1\end{array}\right]} \def\p{\partial} \def\grad#1#2{\frac{\p #1}{\p #2}} \def\qiq{\quad\implies\quad} $Given the vectors $$\eqalign{ h &= \m{h_i\\h_2\\h_3} \quad\qquad x = \m{x_1\\x_2} \\ }$$ the components of the gradient matrix $\LR{G=\LR{\n h}^T=\grad hx}$ are $$\eqalign{ G_{ij} = \grad{h_i}{x_j} \\ }$$ while the component-wise Hessians $\LR{H_i= \n^2h_i}$ are $$\eqalign{ H_{ijk} = \grad{^2h_i}{x_j\,\p x_k} \\\\ }$$
$\sf NB$: The reason for the transpose in the gradient definition is due to the way that the differential is calculated using Nabla notation versus Leibniz notation, i.e. $$ dh \;=\; dx\cdot\n h \;=\; \LR{\grad hx}\cdot dx $$