Partial differential equations (PDEs) play an important role in natural sciences; however, they are far from trivial to solve. This work discusses how high-order time stepping schemes can be implemented in Python and applied to automatically solving PDEs with time-dependent boundary conditions. To this end, the Python libraries FEniCS and Irksome are analyzed. FEniCS provides an automated solution of PDEs employing the Finite Element Method, Irksome implements a way of automating Runge-Kutta time-stepping methods which are commonly used for the time discretization step in solving PDEs. However, they are not compatible. The goal of this work is developing a Python implementation of Runge-Kutta methods which is compatible with FEniCS and based on the implementation provided by Irksome. In general, a main challenge when applying Runge-Kutta methods is their high complexity, especially for higher orders. Therefore, this work implements a solution that automates the calculation of PDEs for different methods by simply specifying their butcher tableau. Furthermore, an example application of solving the elastodynamics equation with the generalized-a method is given.