  We use a recently developed technique, which allows to perform few Lanczos
steps on a given wavefunction even for large system sizes, to investigate the
$t{-}J$ model in the physical parameter region and to check the stability of
the BCS d-wave variational wavefunction. Our statistical Lanczos algorithm,
which extends and improves the one Lanczos step, has been extensively tested on
the small L=18 sites cluster where many Lanczos iterations can be performed
exactly. In this case, at doping $\delta\sim 10%$ the BCS wavefunction
represents a very good initial state to achieve extremely accurate energies and
correlation functions with few Lanczos iterations. For large sizes ($L\le 98$)
the behavior is similar: the low-energy d-wave order parameter $P_d$ is weakly
affected by a couple of Lanczos iterations in the low doping $\delta\sim 10%$
region, whereas the energy is considerably lowered. As a further test of our
calculation we have computed the variance of the Hamiltonian $\Delta E_p=(<
\hat{H}^2 >-< \hat{H} >^2)/L^2$ on the BCS wavefunction with $p=0,1,2$ Lanczos
steps. For large $p$, when the Lanczos algorithm converges to the exact ground
state, the variance vanishes exponentially with increasing $p$. The remarkable
reduction of the variance, observed for $p=1,2$ Lanczos steps even for the
largest lattice size considered, suggests a smooth and rapid convergence to the
exact ground state. These results support the existence of off-diagonal
long-range d-wave superconducting order in the two-dimensional $t{-}J$ model.
