You calculated the gradient wrong (outside valid domain). Check your derivative of the penalty function. The partial derivative of
$$\max(x_1^2+x_2^2-1,0)^2$$
wrt $x_1$ is not
$$4 x_1^3.$$
Rather, it should be
$$4 x_1 (x_1^2+x_2^2-1)$$
(when outside the valid domain).
There are plenty of solutions to this very problem being developed because of the usage of gradient descent methods to optimize neural networks.
Depending on the dimension of your state space $x$ there are smarter or more straight forward ways to do it.
If the dimension of $x$ is low, you can do the brute-force search on the grid surrounding the point you are getting stuck around. How big and dense the grid has to be is dependent on the regularity of the function $f$.
If the dimension is greater than 1 you probably also want to escape the saddle points therefore you should analyse the eigenvalues of a Hessian matrix of $f$ nearby the point you are being stuck.
That being said, if the dimensionality of your problem is low, there are more theoretically sound, higher than order 1 optimization techniques which have some theoretical guarantees on convergence and work better than plain gradient descent.
If the dimensionality of your problem is high, you cannot afford the computation of Hessian due to its quadratic complexity. In this case we want to restrict ourselves to first order optimization techniques. The most popular choices here are SGD, Adam and RAdam. If you are concerned about saddle points you can add noise at each iteration and you can do it smartly as presented in the works of M. Jordan.
As a rule of thumb, if you are looking for good enough solution to your problem and the dimensionality of $x$ is absurdly high, then go and read up on some Deep Learning literature.
If you want any theoretical guarantees of your optimization algorithm and the dimensionality of $x$ is still pretty high then read the Physics literature.
If you have no margin of error and the dimension of $x$ is low then read the Mathematical literature.
The topic is vast and I barely scratched the surface here. I hope you can research the best suited solution for your problem yourself now, because stated as is there is no good answer really.
Best Answer
If a point have derivative zero, the steepest descent converges to that point with a constant sequence, if it starts at that stationary point.
Counterexample required in the comment section:
Consider the function $x \mapsto (\cos((\cos(x))^2_{+})-1)x$, with $t_{+} = t, $ if $t \geq 0$, and $t_{+}=0$, if $t< 0.$ This function is constant in a countable quantity of closed intervals. The quantity of strict maximum is countable and are in between these disjoint intervals in the negative axis. Now, it's not hard to find a negative staring point, say $x_0$, close to a strict local maximum such that in the next step it ends in one of that intervals where the function is constant.