The column routine is not robust, nor is it very efficient. It appears that the existence of multiple solutions for the rx_flash solve block is the source of the problems.
In the column routine, each stage is solved in sequence and then the results for an iteration are compared with those from the previous iteration. For the next iteration, I have tried relaxed direct substitution, Wegstein, and generalized dominant eigenvalue methods. I have also tried different guesses for the rx_flash routine. Solutions can be obtained for some, but not all cases. When solutions were obtained, the computation time was on the order of one hour.
It is possible that the rate based method can avoid the multiple solution problem. It may be awhile before I have something new to post.