Say I have state vectors X1 and X2. X1 contains N1 continuous variables and X2 contains N2 binary (0 or 1) variables. X1 evolves in time according to function f1(X1,X2) and the binary variables in X2 have a probablity per unit time of switching states given by function f2(X1, X2, \dot{X1}).

I’m not sure what keywords I should be searching in the docs to work out if this is solvable using built-in DifferentialEquations.jl features? I’d like to take advantage of variable timestep stiff solvers given by DifferentialEquations.jl as the function f1 is quite costly, but I’m not sure if this setup is too non-trivial for the package, and I’ll just have to implement my own solver?

The mixed jump process + ODE examples are demonstration of this kind of thing with Poisson jump processes instead of binary valued ones. See for example:

What you’d do is where your ODEs and then have callbacks (discrete or continuous) which define the jumping behavior on the discrete values. The growing cell example in the callback portion of the documentation might be helpful: