In [Z. Cai and B. C. Shin, SIAM J. Numer. Anal. 40 (2002), 307-318], we developed the discrete first-order system least squares method for the second-order elliptic boundary value problem by directly approximating $H(div){cap}H(curl)-type$ space based on the Helmholtz decomposition. Under general assumptions, error estimates were established in the $L^2;and;H^1$ norms for the vector and scalar variables, respectively. Such error estimates are optimal with respect to the required regularity of the solution. In this paper, we study solution methods for solving the system of linear equations arising from the discretization of variational formulation which possesses discrete biharmonic term and focus on numerical results including the performances of multigrid preconditioners and the finite element accuracy.