I use a computer algebra to plot a function $f(x,y)$ and notice the existence of a cusp point. See figure below.
Now, I want to determine the exact values $x$ and $y$ of the cusp point. I compute the derivative and solve of $y$ as an expression of $x$, then check the value of $x$ for which the derivative is infinity. Is this approach correct?
Below are the expression of $f(x,y)$, derivative of $f$ and solution of $y$. However, $x = 15$ as the x-coordinate of the cusp point does not make sense.
$$\begin{align} f(x,y) &= -24 x^3+x^2 y^2-16 x^2 y+172 x^2-30 x y^2+264 x y-354 x+y^4-24 y^3+174 y^2-872 y+225\\ \frac{\delta~f(x,y)}{dx} &= -72 x^2+2 x y^2-32 x y+344 x-30 y^2+264 y-354 \end{align}$$
Solving $y$ in $\frac{\delta~f(x,y)}{dx}$ gives $\left\{y\to \frac{-3 \sqrt{4 x^3-72 x^2+189 x+189}+8 x-66}{x-15}\right\},\left\{y\to \frac{3 \sqrt{4 x^3-72 x^2+189 x+189}+8 x-66}{x-15}\right\}$.
As achille hui commented, you need to solve for $(x,y)$ $$\frac{\partial f(x,y)}{\partial x} = \frac{\partial f(x,y)}{\partial y} = 0$$ For simplicity, let us write $$F(x,y)=\frac{\partial f(x,y)}{\partial x}=-72 x^2+2 x y^2-32 x y+344 x-30 y^2+264 y-354$$ $$G(x,y)=\frac{\partial f(x,y)}{\partial y}=2 x^2 y-16 x^2-60 x y+264 x+4 y^3-72 y^2+348 y-872$$ and use Newton-Raphson with $x_0=4$, $y_0=1$ taken from the graph.
Using bilinearization, we have $$0=F(x_0,y_0)+(x-x_0)\frac{\partial F(x,y)}{\partial x}+(y-y_0)\frac{\partial F(x,y)}{\partial y}$$ $$0=G(x_0,y_0)+(x-x_0)\frac{\partial G(x,y)}{\partial x}+(y-y_0)\frac{\partial G(x,y)}{\partial y}$$ where the partial derivatives are computed at $(x_0,y_0)$
This should give as first iterates $$x_1\approx 3.98788 \qquad y_1\approx 1.13939$$ $$x_2\approx 4.00598 \qquad y_2\approx 1.19940$$ $$x_3\approx 4.01485 \qquad y_3\approx 1.22965$$ $$x_4\approx 4.01930 \qquad y_4\approx 1.24482$$ Continuing iterations, you would converge to $$x_\infty\approx 4.02375 \qquad y_\infty\approx 1.26001$$
Edit after achille hui's comment
$F(x,y)$ being a quadratic in $x$ $$F(x,y)=0 \implies x_\pm=\frac{1}{72} \left(y^2-16y+172\pm\sqrt{y^4-32 y^3-1560 y^2+13504 y+4096}\right)$$ Plugging this in $G(x,y)$, did not get anything with $x_-$. However, using $x+$, we end with $$y^5-40 y^4+3880 y^3-63440 y^2+174056 y-447488$$ $${(y^3-24 y^2-780 y+3376) }\,\sqrt{y^4-32 y^3-1560 y^2+13504 y+4096 }=0$$ Then $$\frac{ y^5-40 y^4+3880 y^3-63440 y^2+174056 y-447488} {y^3-24 y^2-780 y+3376}+ \sqrt{y^4-32 y^3-1560 y^2+13504 y+4096 }=0$$ which, after squaring, leads to $$\left(y^3-24 y^2+435 y-512\right)^2\, \, (2 y^2-20 y+113)=0$$ The quadratic factor does not show real roots. So, we are left with$$y^3-24 y^2+435 y-512=0$$ which shows one real root as given by achille hui in his comment. Back to $x_+$, the result !