Introducing a pair-parameter matrix Mittag–Leffler function, we study the uniqueness and Hyers–Ulam stability to a new fractional nonlinear partial integro-differential equation with variable coefficients and a mixed boundary condition using Banach’s contractive principle as well as Babenko’s approach in a Banach space. These investigations have serious applications since uniqueness and stability analysis are essential topics in various research fields. The techniques used also work for different types of differential equations with initial or boundary conditions, as well as integral equations. Moreover, we present a Python code to compute approximate values of our newly established pair-parameter matrix Mittag–Leffler functions, which extend the multivariate Mittag–Leffler function. A few examples are given to show applications of the key results obtained.