In this paper, a numerical algorithm that approximates suitably the solution of certain non-autonomous partial differential systems of equations is constructed. Firstly, conditions to be satisfied for the data of the problem to ensure the existence of a unique solution is established and, after using a matrix separation of variables technique, an exact series solution is constructed. Each term in the series solution must solve a linear system of non-autonomous differential equations with oscillatory solutions. In order to approximate numerically the formal series solution, schemes based on the Magnus expansions are considered, and we show that commutator-free Magnus integrators are of great interest for solving this problem. We show that this family of methods admit a recursive relation which provides numerical approximations for all terms in the series solution. This iterative recursion has similar error growth as the exact solution, and this ensures the convergence of the numerical solution. These exponential integrators usually often provide good results, generally better than other standard numerical methods, both from the standpoint of preserving qualitative properties of the problem and the computational cost. Finally, numerical experiments are considered in order to test the behavior of the constructed algorithm.