This turned out to be a great solution, thank you!
For those interested, I ended up doing this to get the eigenvalue with largest real part.:
using ForwardDiff, NonlinearSolve
function stability(inits, parameters) #initial conditions and parameters
NLprob = NonlinearProblem(Alpha_Temp_B, inits, parameters) # find roots of model Alpha_Temp_B(u,p)
NLsol = solve(NLprob, NewtonRaphson()) # find the steady-state
J = ForwardDiff.jacobian(u -> Alpha_Temp_B(u, parameters), NLsol.u) # get Jacobian
max_real = maximum(real.(eigen(J).values)) #get eigenvalue with maximum real part
end
I had previously found convenient code to do this here: Incorrect bifurcation diagram bifurcationkit.jl - #4 by dawbarton