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, σ) $ 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, σ) $ 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, σ) $, 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.