Consider the system NDSolve[{y'[x]==g[y[x],x],y[0]==q},y[x],{x,0,xf}]; g[x,y[x]]=Module[{gs,res,g0=0},gs=FindRoot[h[x,y[x]]==0,{g.g0}];res=g/.gs[[1]];Return[res]
This always gives error messages of the form FindRoot::nlnum: "The function value {8.51808*10^10-1.\ (0.0333333\x+0.333333\y[x])^(2/3)\ (0.05\x+y[x])^(1/3)} is not a list of numbers with dimensions {1} at {p2$97812} = {0.05`}." where all the above numbers are coming from a fairly complicated expression for h
It seems clear to me that NDSolve regards g[x,y[x]] as an analytical function rather than the numerical result of a calculation. It is probably doing this in order to calculate the derivative of g with respect to x.
If this is what is going on, is there some way to force NDSolve to use the chord method to calculate the derivative of g with respect to x? Only FindRoot seems to be able to get a numerical value of g from the expression for h.
