We present an O(n) complexity and implicit algorithm for the two-dimensional solution of the Stream Power Incision Model (SPIM) enriched by a discharge threshold term and taking into account variability in rainfall and thus discharge. The algorithm is based on the formulation developed by Deal et al (2018) and the generalization of the FastScape algorithm \cite{BraunWillett2013}where the slope is approximated by first-order accurate finite difference. We consider a variety of discharge thresholds that vary in their dependence on channel slope. The algorithm requires finding the root of a non-linear equation using a Newton-Raphson iterative scheme. We show that the convergence of this scheme is unconditional, except for a narrow range of model parameters where the threshold increases with the slope and for low discharge variability. We also show that the rate of convergence of the iterative scheme is directly proportional to the slope exponent n in the SPIM. We compare the algorithm to analytical solutions and to numerical solutions obtained using a higher-order finite difference scheme. We show that the accuracy of the FastScape algorithm and its generalization presented here is comparable to other schemes for values of n>1. We also confirm that the FastScape algorithm and its generalization to variable discharge+threshold conditions does not need to satisfy the CFL condition and provides an accurate solution for both small and very long time steps. We finally use the new algorithm to quantify how the existence of an erosional threshold strongly affects the length of the post-orogenic decay of mountain belts.