  We integrate numerically the Kardar-Parisi-Zhang (KPZ) equation in 1+1 and
2+1 dimensions using an Euler discretization scheme and the replacement of
${(\nabla h)}^2$ by exponentially decreasing functions of that quantity to
suppress instabilities. When applied to the equation in 1+1 dimensions, the
method of instability control provides values of scaling amplitudes consistent
with exactly known results, in contrast to the deviations generated by the
original scheme. In 2+1 dimensions, we spanned a range of the model parameters
where transients with Edwards-Wilkinson or random growth are not bserved, in
box sizes $8\leq L\leq 128$. We obtain roughness exponent $0.37\leq \alpha\leq
0.40$ and steady state height distributions with skewness $S= 0.25\pm 0.01$ and
kurtosis $Q= 0.15\pm 0.1$. These estimates are obtained after extrapolations to
the large $L$ limit, which is necessary due to significant finite-size effects
in the estimates of effective exponents and height distributions. On the other
hand, the steady state roughness distributions show weak scaling corrections
and evidence of stretched exponentials tails. These results confirm previous
estimates from lattice models, showing their reliability as representatives of
the KPZ class.
