Jumping in on your exchange: using the solver does not prevent Stan from accurately calculating gradients. There are pathological cases where the gradients will not be well defined, for instance when your equation is degenerate, or when f(z, \theta) is discontinuous with respect to either the unknown z or the parameters \theta. In such a case, your posterior is discontinuous and you pick up very distinct symptomes (which I believe are not observed here).
I ran into such a problem with a model which embedded a KKT problem (an algebraic equation with an equality and an inequality constraint). See Non-smooth posterior and KKT problem.