For example I have a function:
$$f(x_1,x_2) = \begin{bmatrix} x_1x_2^2 + x_1^3x_2\\x_1^2x_2 + x_1 + x_2^3\\\end{bmatrix}$$
Is it possible to find a gradient of this function? Because knowing the interpretation of gradient(direction of the steepest slope) i think it is not possible. We can do it, if function has codomain $R^1$
Am I right? If wrong, explain me please.
Your function $f(x_1,x_2)$ is of the form $F(x_1,x_2) = f(g(x_1,x_2),h(x_1,x_2))$ I think. Where $g$ and $h$ are the first and second entry respectively. So the gradient that you are looking for is of the form:
$J_f(x_1,x_2) = \begin{bmatrix} \frac{\partial{g(x_1,x_2)}}{\partial{x_1}} & \frac{\partial{g(x_1,x_2)}}{\partial{x_2}} \\ \frac{\partial{h(x_1,x_2)}}{\partial{x_1}} & \frac{\partial{h(x_1,x_2)}}{\partial{x_2}} \\\end{bmatrix}$ I think that in order to find the directional derivative in the direction of a unit vector $\bar{u} = (u_1,u_2)$, you just need to multiply the matrix and the vector: $J_f(x_1,x_2)\bar{u}^T$