  We integrate for the first time the hydrodynamic
Hall-Vinen-Bekarevich-Khalatnikov equations of motion of a $^{1}S_{0}$-paired
neutron superfluid in a rotating spherical shell, using a pseudospectral
collocation algorithm coupled with a time-split fractional scheme. Numerical
instabilities are smoothed by spectral filtering. Three numerical experiments
are conducted, with the following results. (i) When the inner and outer spheres
are put into steady differential rotation, the viscous torque exerted on the
spheres oscillates quasiperiodically and persistently (after an initial
transient). The fractional oscillation amplitude ($\sim 10^{-2}$) increases
with the angular shear and decreases with the gap width. (ii) When the outer
sphere is accelerated impulsively after an interval of steady differential
rotation, the torque increases suddenly, relaxes exponentially, then oscillates
persistently as in (i). The relaxation time-scale is determined principally by
the angular velocity jump, whereas the oscillation amplitude is determined
principally by the gap width. (iii) When the mutual friction force changes
suddenly from Hall-Vinen to Gorter-Mellink form, as happens when a rectilinear
array of quantized Feynman-Onsager vortices is destabilized by a counterflow to
form a reconnecting vortex tangle, the relaxation time-scale is reduced by a
factor of $\sim 3$ compared to (ii), and the system reaches a stationary state
where the torque oscillates with fractional amplitude $\sim 10^{-3}$ about a
constant mean value. Preliminary scalings are computed for observable
quantities like angular velocity and acceleration as functions of Reynolds
number, angular shear, and gap width. The results are applied to the timing
irregularities (e.g., glitches and timing noise) observed in radio pulsars.
