The discrete nonlinear Schrödinger (DNLS) equation is a nonintegrable dynamical lattice model with complex behavior and a growing list of applications in quantum field theory, condensed matter physics, and quantum optics. In this thesis, classical soliton solutions of the DNLS are studied numerically in a version of the model with cubic and quintic nonlinear interactions. Starting with initial soliton profiles that are obtained from a semi-analytical variational approximation, or by iteration on the homoclinic orbits of a two-dimensional map, numerically exact initial profiles are generated by using a multi-dimensional Newton iteration algorithm. The stability of these solutions is investigated in appropriate ranges of the parameter space of the model. From these stationary exact soliton solutions, initial profiles for moving soliton solutions are constructed by deformation with a complex phase function. Several different numerical algorithms have been proposed to integrate time-dependent equations like the DNLS: three of the most stable and efficient methods are used to obtain moving soliton solutions, and their accuracies and efficiencies are compared. In particular, an adaptive symmetrized split-step Fourier method that was very recently proposed to solve continuum soliton equations is modified and applied to this DNLS problem. Finally, collisions of two moving solitons are studied using the same integration algorithms, with results that illustrate a wide range of complex behavior that can result from the characteristic discreteness and nonintegrability of the DNLS.