The Trilinos library LOCA (http://www.cs.sandia.gov/LOCA/ ) allows computing branches of steady states of large-scale dynamical systems like (discretized) nonlinear PDEs. The core algorithms typically are (pseudo-)arclength continuation, Newton–Krylov methods and (sparse) eigenvalue solvers. While LOCA includes some basic techniques for computing bifurcation points and switching branches, the exploration of a complete bifurcation diagram still takes a lot of programming effort and manual interference. On the other hand, recent developments in algorithms for fully automatic exploration are condensed in PyNCT (https://pypi.org/project/PyNCT/ ). The scope of this algorithmically versatile software is, however, limited to relatively small (e.g. 2D) problems because it relies on linear algebra from Python libraries like NumPy. Furthermore, PyNCT currently does not support problems with a non-Hermitian Jacobian matrix, which rules out interesting applications in chemistry and fluid dynamics. In this paper we aim to combine the best of both worlds: a high-level implementation of algorithms in PyNCT with parallel models and linear algebra implemented in Trilinos. PyNCT is extended to non-symmetric systems and its complete backend is replaced by the PHIST library (https://bitbucket.org/essex/phist ), which allows us to use the same underlying HPC libraries as LOCA does. We then apply the new code to a reaction-diffusion model to demonstrate its potential of enabling fully automatic bifurcation analysis on parallel computers.