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
$$\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)$$
