The coherent state variational algorithm provides a method for solving the large N limit of non-abelian gauge theories. An implementation of this algorithm, capable of minimizing the large N effective action and computing meson and glueball spectra, has recently been completed. Hamiltonian or euclidian formulations of lattice gauge theories, in any dimension, may be studied. Bose or fermi fundamental representation matter fields may be included. This paper discusses the design and testing of this implementation. The method involves explicit manipulation of expectation values of physical operators and may be applied directly in infinite volume. The error introduced by the truncation of the set of physical observables (necessary to obtain a finite procedure) is studied by applying the algorithm to a variety of exactly soluble model theories. These include φ4 scalar field theories, (ψψ)2 fermion theories, 2-dimensional euclidean pure gauge theory, and 1+1 dimensional QCD. Modest size calculations are shown to yeild accurate results, even in theories possessing asymptotic freedom, spontaneous symmetry breaking, or large N phase transitions.