For multiple-input multiple-output (MIMO) spatial-multiplexing transmission, zero-forcing detection (ZF) is appealing because of its low complexity. Our recent MIMO ZF performance analysis for Rician--Rayleigh fading, which is relevant in heterogeneous networks, has yielded for the ZF outage probability and ergodic capacity infinite-series expressions. Because they arose from expanding the confluent hypergeometric function $ {_1! F_1} (cdot, cdot, sigma) $ around 0, they do not converge numerically at realistically-high Rician $ K $-factor values. Therefore, herein, we seek to take advantage of the fact that $ {_1! F_1} (cdot, cdot, sigma) $ satisfies a differential equation, i.e., it is a textit{holonomic} function. Holonomic functions can be computed by the textit{holonomic gradient method} (HGM), i.e., by numerically solving the satisfied differential equation. Thus, we first reveal that the moment generating function (m.g.f.) and probability density function (p.d.f.) of the ZF signal-to-noise ratio (SNR) are holonomic. Then, from the differential equation for $ {_1! F_1} (cdot, cdot, sigma) $, we deduce those satisfied by the SNR m.g.f. and p.d.f., and demonstrate that the HGM helps compute the p.d.f. accurately at practically-relevant values of $ K $. Finally, numerical integration of the SNR p.d.f. produced by HGM yields accurate ZF outage probability and ergodic capacity results.