<p>Recently, the algebraic method of Mikusiński’s operational calculus was extended for the first time to fractional PDEs on orthant domains. We now investigate how the method can be adapted to solve fractional PDEs involving Caputo derivatives, with more natural boundary conditions than the previous work with Riemann–Liouville derivatives. We create an algebraic framework for interpreting multi-dimensional integro-differential operators and solve some PDEs in two dimensions, with the solutions being written in terms of some newly defined multivariate Mittag-Leffler functions.</p>