OK, here's a quick rundown on how I would do this; hopefully our OP Mirjan Pecenko can check his own work against what I do here.
With
$\dot x = xy + 1, \tag 1$
$\dot y = x + xy, \tag 2$
we have equilibria where
$\dot x = xy + 1 = \dot y = xy + x = 0; \tag 3$
thus,
$xy = -x \Longrightarrow -x + 1 = 0 \Longrightarrow x = 1, \tag 4$
and then
$xy + 1 = 0 \Longrightarrow xy = -1 \Longrightarrow y = (1)y = -1; \tag 4$
so it appears the only critical point is at $(1, -1)$. The Jacobean of the vector field
$\vec V(x, y) = \begin{pmatrix} xy + 1 \\ x + xy \end{pmatrix} \tag 5$
at $(1, -1)$ is
$J_V(1, -1) = \begin{bmatrix} (\partial(xy + 1)/\partial x & (\partial(xy + 1)/\partial y \\ (\partial(xy + x)/\partial x & (\partial(xy + x)/\partial y \end{bmatrix}_{(1, -1)}$
$= \begin{bmatrix} y & x \\ y + 1 & x \end{bmatrix}_{(1, -1)} = \begin{bmatrix} -1 & 1 \\ 0 & 1 \end{bmatrix}; \tag 6$
it is now obvious that the eigenvalues of $J_V(1, -1)$ are $\pm 1$; therefore this point is a saddle; the eigenvectors at $(1, -1)$ are $(1, 0)^T$ for $-1$ and $(1/2, 1)^T$ for $1$; it is now easy to sketch a phase portrait for this system, a task I leave to my readers.
It should be noted that when sketching a phase portrait, it is often helpful to find those curves in $\Bbb R^2$ where $\dot x = 0$ and/or $\dot y = 0$. These curves are useful guides to finding the geometry of the solutions, since they show us where the tangent lines to the integral curves or vertical or horizontal, respectively. When combined with the shapes of the solutions near the critical point provided by the above analysis, we can get a pretty good idea of how the flow will appear. As with any hand-done graphical analysis, care must be taken to ensure that we draw accurately enough to capture only correct features of the trajectories.
Note Added in Edit, Thursday 14 June 2018 12:35 PM PST: It appears that the notion of isoclines, which proves to be most convenient as a guide to sketching phase portraits, may be generalized in a way which allows the extraction of more information about the integral curves and/or overall shape of the flow of a given vector field. In this problem, isoclines are only exploited in a rather peripheral way, since they are merely mentioned as a sort of after-thought in the comments. Nevertheless, they may be used much more extensively. Indeed, rather than restricting the use of the isocline method to determining the curves on which $\dot x = 0$ and/or $\dot y = 0$, we may if we so choose deploy it in an attempt to find just where $\dot x, \dot y$ take on any of their possible values. One technique which can help effect this aim is to use the gradient of the component functions, in this case $xy + 1$ and $xy + x$, to guide us towards regions of greater or lesser component magnitude; for example, since
$\dot x = xy + 1, \tag 7$
we have
$\nabla (\dot x) = (y, x)^T, \tag 8$
which indicates the direction in which $\dot x$ increases, so that we may, for example, find the directions in which points on the $\dot x = 0$ isocline must be moved to make $\dot x$ larger. By systematic application of such methodology, quite detailed phase portraits may be had. Unfortunately, at present I lack both the graphics tools and the time to illustrate what I am saying via pictorial means. End of Note.
