An algorithm is developed for fully nonlinear three-dimensional (3D) simulation of a difference-frequency acoustic beam resulting from the interaction of two high-intensity pump waves. Simulations are performed in the frequency domain based on the Khokhlov-Zabolotskaya-Kuznetsov equation. A spectrum filtering method is used to enable accurate solutions for the difference-frequency fields in strongly nonlinear beams and with a high downshift frequency ratio using only dozens of spectral components retained in the algorithm. As an example, the dual-frequency operation of an underwater multi-element ellipsoidal array is considered, and numerical solutions describing parametric interactions in the array field are analyzed. It is shown that difference-frequency beams are more symmetric in transverse directions compared with the pump beams. The most efficient parametric generation of difference-frequency beams corresponded to close and beyond shock-forming conditions. Axial pressure amplitude of the difference frequency was shown to grow first quadratically with the source pressure following the quasi-linear solution and then linearly once shocks start to develop. The percentage of the total power converted to the difference frequency from pump waves increased at high power outputs without saturation. Up to twofold increase in directivity angles of difference-frequency beams under shock-forming conditions was observed compared with quasi-linear conditions.