Most methods (that I know of) of numerically approximating the solution of ODEs are "general linear methods". For this type of method, the so-called 'linear stability' is examined by applying the method to the linear, constant-coefficient, (complex) scalar ODE

$\dot{y} = \lambda y $.

Because the same procedure is used to generate a new approximation at each time step, the methods result in linear, homogeneous, constant-coefficient difference equations for the approximate solution values at each time step. As a result, linear stability analysis for these methods amounts to analyzing the solutions of constant coefficient, homogeneous linear difference equations.

If you vary the method randomly at each time step, the coefficients of the difference equation are no longer constant. For example, second order explicit linear multistep methods can be written as a one-parameter family. If you choose this parameter randomly at each time step, the difference equation looks like this

$y_{n+1} + F(a(n))y_n + G(a(n))y_{n-1} = 0$

where $a(n)$ is a random variable (the parameter in the family of methods) and the functions $F$ and $G$ are known. My main question is whether there is any theory giving conditions on $F$, $G$, and the distribution of $a$ such that solutions remain bounded in the limit of large $n$. It would be nice if the theory generalized to higher order difference equations too.

If you were to randomly decide to use second-order multistep versus second-order Runge-Kutta methods, for example, then $F$ and $G$ would also depend on $n$. A theory to handle that case would be welcome too.

Since the numerical analyst is free to choose the distribution of $a$, and to some extent also $F$ and $G$, I'm wondering if it might be possible to design 'random' methods that have better linear stability properties than the usual ones which repeat the same process over and over. I'm posting here because I know almost nothing about stochastic/random processes.