  The purpose of this work is to introduce a new idea of how to avoid the
factorization of large matrices during the solution of stiff systems of ODEs.
Starting from the general form of an explicit linear multistep method we
suggest to adaptively choose its coefficients on each integration step in order
to minimize the norm of the residual of an implicit BDF formula. Thereby we
reduce the number of unknowns on each step from $n$ to $O(1)$, where $n$ is the
dimension of the ODE system. We call this type of methods Minimal Residual
Multistep (MRMS) methods. In the case of linear non-autonomous problem, besides
the evaluations of the right-hand side of ODE, the resulting numerical scheme
additionally requires one solution of a linear least-squares problem with a
thin matrix per step. We show that the order of the method and its
zero-stability properties coincide with those of the used underlying BDF
formula. For the simplest analog of the implicit Euler method the properties of
linear stability are investigated. Though the classical absolute stability
analysis is not fully relevant to the MRMS methods, it is shown that this
one-step method is applicable in stiff case. In the numerical experiment
section we consider the fixed-step integration of a two-dimensional
non-autonomous heat equation using the MRMS methods and their classical BDF
counterparts. The starting values are taken from a preset slowly-varying exact
solution. The comparison showed that both methods give similar numerical
solutions, but in the case of large systems the MRMS methods are faster, and
their advantage considerably increases with the growth of dimension. Python
code with the experimantal code can be downloaded from the GitHub repository
https://github.com/bfaleichik/mrms.
