  In the modelling of blazars and other radiating plasma flows, correctly
evolving in time a series of kinetic equation coupled for different particle
species is a computationally challenging problem. In this paper we introduce
the code Katu, capable of simultaneously evolving a large number of particle
species including photons, leptons and hadrons and their interactions. For each
particle, the numerical expressions for the kinetic equations are optimized for
efficiency and the code is written to easily allow for extensions and
modifications. Additionally, we tie Katu to two pieces of statistical software,
emcee and pymultinest, which makes a suite ready to apply to blazars and
extract relevant statistical information from their electromagnetic (and
neutrino, if applicable) flux. As an example, we apply Katu to Mrk 421 and
compare a leptonic and a leptohadronic model. In this first ever direct
Bayesian comparison between full kinetic simulations of these models, we find
substantial evidence (Bayes factor $\sim7$) for the leptohadronic model purely
from the electromagnetic spectrum. The code is available online at
https://github.com/hveerten/katu under the LGPL license.
