The zero-lag state covariance matrix of a mechanical system subject to white-noise excitation is computed by algebraic techniques. This matrix is of interest because it is the solution to the familiar mean-square response computations. The problem is to determine the solution of the matrix equation AP+PA′= −G for P given A, an n×n nonsingular (structural system) matrix, and G, a symmetric, positive, semidefinite, input covariance matrix. P is the covariance matrix in question. This problem is the stochastic dual of the direct method of Lyapounov in the study of the stability of linear systems. The equation is solved in two ways. The direct solution involves inversion of a matrix whose order increases as the square of the number of degrees-of-freedom (df). Another solution is accomplished by a decomposition in which the input excitation is separated into stochastically independent components. In that case, for each input, G has rank 1 and the solution is readily computable upon transformation to a convenient form by a technique introduced by Kryloff in 1937. In the vernacular of contemporary control theory, the problem solution employs canonical transformations and the concept of a controllability matrix introduced and exploited by Kalman in 1960. The second solution of the equation requires inversion of a matrix whose order increases linearly with the number of degrees-of-freedom. The algebraic techniques discussed here are simpler to use than the “classical,” complex, variable-residue-theory techniques. More important, and in contrast with the classical techniques, the procedures described here are applicable to systems with an arbitrarily large number of degrees of freedom and may also be used to compute mean-square responses for systems with a vector of inputs. Explicit computational formulas for mean-square responses computed by the direct method are given for the general 1-, 2-, and 3-df systems. Numerical and graphical examples illustrate the procedure. Computational results, procedures, and formulas for the alternate “indirect” method are also worked for the 1-, 2-, and 3-df systems. Machine computations on a 2-df system illustrate the details for the practical computational use of the procedure.