***
> _Or how else is such a non-linear relationship without closed-form
> solution approached?_


For maximum likelihood we want the density of $y$ conditional on $x$ (because this is our sample), $f_{y|x}(y|x)$.

We start with making an assumption on the density of $\varepsilon$ conditional on $x$, $f_{\varepsilon|x}(\varepsilon|x)$.

When the relation is not implicit, say $y = ax + \varepsilon$ the change-of-variables method has a Jacobian determinant equal to unity so we simply have

$$f_{y|x}(y|x) = f_{\varepsilon|x}(y-ax|x)$$

and we can proceed as usual.
But when the relation is implicit we have

$$\varepsilon = y-h(y,\theta)-x \implies \frac{\partial \varepsilon}{\partial y} = 1- \frac{\partial h(y,\theta)}{\partial y}$$

so here, the observation density will be

$$f_{y|x}(y|x) = \left|1- \frac{\partial h(y,\theta)}{\partial y}\right| \cdot f_{\varepsilon|x}(y-h(y,\theta)-x|x)$$