A direct method for the computation of polynomial conservation laws of polynomial systems of nonlinear partial differential equations (PDEs) in multi-dimensions is presented. The method avoids advanced differential-geometric tools. Instead, it is solely based on calculus, variational calculus, and linear algebra.Densities are constructed as linear combinations of scaling homogeneous terms with undetermined coefficients. The variational derivative (Euler operator) is used to compute the undetermined coefficients. The homotopy operator is used to compute the fluxes.The method is illustrated with nonlinear PDEs describing wave phenomena in fluid dynamics, plasma physics, and quantum physics. For PDEs with parameters, the method determines the conditions on the parameters so that a sequence of conserved densities might exist. The existence of a large number of conservation laws is a predictor for complete integrability. The method is algorithmic, applicable to a variety of PDEs, and can be implemented in computer algebra systems such as Mathematica, Maple, and REDUCE.