For those of you familiar with the method, it is known that one must solve the equation:
y(i+1) = y(i) + h*F( X(i+1), Y(i+1) )
However, F is usually not linear, and the resulting equation usually has many different solutions for y(i+1). Which solution are we looking for, and what would one do, for instance, to the Newton-Raphson method so it finds the CORRECT zero? Any help is appreciated.
Usually you chose the solution that is closest to y(i). The correct solution should satisfy
y(i+1)=y(i)+h*f(x(i),y(i)) + O(h²)
However, for stiff problems the constant in O(h²) can be very large so that this relation is not as helpful as it seems.
If L is the Lipschitz constant of f and q=L·h < 1 then the iteration
k(j+1) = f(x(i+1), y(i)+h*k(j))
converges geometrically with error proportional to q^j/(1-q) or smaller.
This is violated for stiff problems, where you use the Jacobian J of f in (x(i+1),y(i)) to iterate
k(j+1) = k(j) - (I-h*J)^(-1)*( k(j)-f( x(i+1), y(i)+h*k(j)) )
possibly approximating (I-h*J)^(-1) as I+h*J, see Rosenbrock method. Starting with k(0)=0 or k(0)=f(x(i),y(i)) should almost always rapidly converge to the correct solution.