We need to decide what external library we need to use (if any) - mainly for the solver. In @Kazu1886's original code, he used a GMRES solver provided by the GSL library.
Do we have to use an external solver? Can the equations be rewritten and solved without an external solver?
Is GMRES the best solver for this model?
What is the best library to use? I found this Fortran interface to GSL, so I guess we can try this? A more powerful library that includes nonlinear iterative solvers seems to be PETSc.
We need to decide what external library we need to use (if any) - mainly for the solver. In @Kazu1886's original code, he used a GMRES solver provided by the GSL library.