A method of constructing a numerical theory of comet motion covering large intervals of time is presented. It involves the determination of individual values of the constants Α₁, Α₂, Α₃ (radial, transverse, and normal components of non-gravitational acceleration) and the displacement of the photo center in each apparition. In difficult cases, such as close approaches to the major planets, sudden increase of brightness, or considerable asymmetric outgassing about perihelion, when Marsden's model of non-gravitational acceleration does not allow one to present observations with sufficient accuracy, the instantaneous variation of velocity in the cometary motion is as-sumed. Based on this methodology, a unified numerical theory of the motion of comet Kopff during the time interval 1906-2002 was developed. It links 16 cometary apparitions with an RMS residual σ = 1.40".