Hi, I'm new to posting questions on forums, so I apologise if the problem is poorly described. My problem is solving a simulation of the state of a mineral processing froth flotation plant. In the form x@i+1 = f(x@i), f represents the flotation plant. f is a computationally intensive solution of systems of simultaneous linear equations (material balancing) and evaluation of the state of each unit in the process. x contains ~20 elements, depending on user configuration of the process. For some (user defined) configurations of the function f, this system iterates rapidly to convergence by back-substitution ie. Picard iteration. For the last couple of years I have been using Krasnoselskij iteration (EMA filter) and the system converges in most, but not all situations. Also, convergence is slow (200+ iterations) for some configurations. Since f is computationally expensive, on the order of 10ms, calculation of a good x@i+1 is crucial. I have also investigated Mann iterations. I have implemeted Direct Inversion in the Iterative Subspace (DIIS) but still have concerns about the uniqueness of solutions from DIIS which I am investigating. A single calculation of the Jacobian would take longer than most simulations take to converge by back-substitution. Is there a recommended algorithm for accelerating convergence and increasing the stability of a multidimensional fixed point iteration problem without derivatives?