A covariant shell model of three Dirac fermions in the center-of-mass frame is analytically solved for a particular spin-dependent harmonic-oscillator potential. Three equal-mass quarks are all assumed to be in states with identical quantum numbers ofjπ=1/2+ for the shell-model ground state. The sixtyfour-component composite wave function for three Dirac fermions is reduced to a set of four coupled differential equations for the four independent composite radial components. Hyperspherical coordinates are used. The two-body potential used vanish in the small component of the composite wave function. The various components of the composite wave function are found to have aρl exp(−κρ2) hyperradial dependence, where ρ is the hyperradius, andl is the sum of the orbital angular momentum acting in the various components of the composite wave function. For an eigenergy equal to the proteon rest-mass energy, the root-mean-square radius calculated is about 0.84 fermi, about 80% of that calculated in an independent-particle harmonic-oscillator model with the same rest-mass energy for the system. A sixth-order polynomial equation is found for the rest-mass energy of the system in terms of the potential-strength parameter, and the quark mass, For a quark mass of zero, a quadratic equation inE3 is obtained. One solution agrees with earlier work, and the other can be rejected on physical grounds. The various components of the composite three-body wave function are dominated by thel=0s-state component as the quark mass approaches one third of the system rest-mass energy. For zero quark mass thes-wave component is about 29% of the normalization, thep wave about 44%, thed wave about 22%, and thef wave 4% of the normalization.