First order Non-Linear ODE with no explicit form of the derivative I have the following non-linear first order ODE
$$
[a+b(1-e^{-m\frac{dy}{dx}})]\frac{dy}{dx}=f(x)
$$
to be integrated over the range $x=x_0$ to $x=x_1$
The ODE is of the form
$$
y'=f(x,y,y')
$$
which does not fit into the framework of the Runge-Kutta algorithm.
The approach that I had in mind was to subdivide the interval $[x_0,x_1]$ into subintervals and in each subinterval, set $\phi=\frac{dy}{dx}$ and solve the non-linear equation
$$
[a+b(1-e^{-m\phi})]\phi-f(x)=0
$$
at the start and end-points of the sub-interval and at two points in the middle, all equally spaced. With $\phi$ known at these 4 points, I could fit a cubic polynomial to the derivative as
$$
\frac{dy}{dx}=A+Bx+cx^2+Dx^3
$$
and obtain the value at the end point of the sub-interval as
$$
y_{i}=y_{i-1}+A_i(x_i-x_{i-1})+\frac{B_i}{2}(x_i^2-x_{i-1})
+\frac{C_i}{3}(x_i^3-x_{i-1}^3)+\frac{D_i}{4}(x_i^4-x_{i-1}^4)
$$
The question that I have is whether there is a more sophisticated approach that I should consider, possibly even some variant of traditional RK methods.
Thanks
 A: Essentially, you are given $dy/dx$ (though implicitly, but you can calculate it wherever you need it), so it's just a question of numerical integration. There's no need to find fitting polynomials, you can just plug those function values into known quadrature rules exact for cubic polynomials, like Simpson't 3/8 rule. If you need a list of values of $y(x)$ for some equispaced $x$, that may be the most practical approach. Details depend on the regularity of your functions (including $f(x)$) and the values of your parameters (if $m$ is large or very large, there may be problems).
Concerning the "more sophisticated approach": depending on the shape of your RHS $f(x),$ there may be an analytical solution in a parametric way, $x=x(\phi), y=y(\phi).$ Just a very simple case meant only as an illustration:
Let's consider the equation $$\left(2-e^{-\frac{dy}{dx}}\right)\frac{dy}{dx}=x,$$ with initial condition $y(0)=0.$ $x$ is a function of $\phi=\frac{dy}{dx},$ $$x(\phi)=\left(2-e^{-\phi}\right)\phi,$$ so we can find $y$ as a function of $\phi,$ too: $$\frac{dy}{d\phi}=\phi\frac{dx}{d\phi}=2\phi-(\phi+\phi^2)\,e^{-\phi},$$ i.e.
$$y(\phi)=\phi^2+(\phi^2+3\phi+3)\,e^{-\phi}-3,$$ considering the initial condition (as $x=0$ corresponds to $\phi=0$). That's sufficient for a parametric plot of the solution:

It's not very impressive, I'll admit, just a proof of concept.
