I am trying to implement a stiff ODE solver that uses step doubling with the backward Euler method but I want to make a modification and I am unsure how to incorporate it.I am trying to solve $y'=f(y), y(a)=d$ and I am using an approximation to $f$ such that $y'\approx F(y),\;t\in (T,T+dt), y(T)=Y$ where$$F(y)=f(y)+Df(Y)(y(t) - Y).$$I understand how to implement the backward Euler method when we simply have $y'=f(t,y)$ but I am having trouble writing a program to solve the problem here because $F(y)$ is dependent on $y(t)$, which is obviously unknown so the usual formulation of $\Delta x=Df^{-1}(x)f(x)$ doesn't seem to work. I was wondering how I could go about writing a method for this? 
Slight Modification to Backward Euler Stiff ODE Solver
ode;nonlinear equations;implicit methods
Let's assume you have some differential equation that can be written as the following:$$ \frac{dy}{dt} = f(t,y)$$for some initial condition $y(a) = d$. We then want to integrate this quantity using Backward Euler, with the fundamental formula being derived in the following way:$$ \begin{align}\frac{dy}{dt} &= f(t_{n+1},y_{n+1})\\\frac{y_{n+1} - y_{n}}{\Delta t} &= f(t_{n+1},y_{n+1})\\y_{n+1} &= y_{n} + \Delta t f(t_{n+1},y_{n+1})\end{align}$$Now we might want to make things a little easier than solving nonlinear equations for $y_{n+1}$, if we're willing to add in a little bit of error, by linearizing $f(t_{n+1},y_{n+1})$ around $y_{n}$. We can show the linearization for $f(t_{n+1},y_{n+1})$ about $y_{n}$ is:$$ f(t_{n+1},y_{n+1}) = f(t_{n+1},y_{n} + (y_{n+1} - y_{n})) \approx f(t_{n+1},y_{n}) + \left.\frac{\partial f}{\partial y}\right|_{(t_{n+1},y_{n})}\left(y_{n+1} - y_{n}\right)$$If we substitute this into the Backward Euler equation, and set $A_n = \left.\frac{\partial f}{\partial y}\right|_{(t_{n+1},y_{n})}$, we can work out the following:$$ \begin{align}y_{n+1} &= y_{n} + \Delta t f(t_{n+1},y_{n}) + \Delta tA_n\left(y_{n+1} - y_{n}\right) \\(I + \Delta t A_n)y_{n+1} &= (I - \Delta t A_n)y_{n} + \Delta t f(t_{n+1},y_{n}) \\y_{n+1} &= y_{n} + \Delta t (I - \Delta t A_n)^{-1}f(t_{n+1},y_{n})\end{align}$$The final formula above should be all you need, assuming you can compute $A_n$.