We have developed a code for modeling 3-dimensional kinematical data (i.e., both radial velocity and proper motions) of individual tracers in stellar systems. The code relies on the Schwarzschild method of orbit superposition to find the distribution function (a function of energy, angular momentum, and third integral of motion) that best reproduces all the observations (the positions and velocities of the tracers as well as the overall light distribution) in a given potential. Previous implementations of this method have generally focused on fitting integrated light measurements, instead of data for individual tracers. However, the latter is the only data that are generally available in a large class of problems (e.g., kinematics of tracers in nearby globular clusters, kinematics of planetary nebulae or globular clusters in galaxy halos, etc.). We have performed extensive tests of the fitting process by making use of data generated by a well-tested N-body simulation of an axisymmetric galaxy, quantifying how well the code can recover the known input distribution function and potential.