  We investigate numerical solution of $Q^2$ evolution equations for structure
functions in the nucleon and in nuclei.
(Dokshitzer-Gribov-Lipatov-)Altarelli-Parisi and Mueller-Qiu evolution
equations are solved in a brute-force method. Spin-independent
flavor-nonsinglet and singlet equations with next-to-leading-order $\alpha_s$
corrections are studied. Dividing the variables $x$ and $Q^2$ into small steps,
we simply solve the integrodifferential equations. Numerical results indicate
that accuracy is better than 2\% in the region $10^{-4}<x<0.8$ if more than
two-hundred $Q^2$ steps and more than one-thousand $x$ steps are taken. The
numerical solution is discussed in detail, and evolution results are compared
with $Q^2$ dependent data in CDHSW, SLAC, BCDMS, EMC, NMC, Fermilab-E665, ZEUS,
and H1 experiments. We provide a FORTRAN program for Q$^2$ evolution (and
``devolution'') of nonsinglet-quark, singlet-quark, $q_i+\bar q_i$, and gluon
distributions (and corresponding structure functions) in the nucleon and in
nuclei. This is a very useful program for studying spin-independent structure
functions.
