For multiple-input–multiple-output (MIMO) spatial-multiplexing transmission, zero-forcing (ZF) detection 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 <bold>holonomic</bold> function. Holonomic functions can be computed by the <bold>holonomic gradient method</bold> (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.