# Get Started

using FindSteadyStates
using DifferentialEquations

## Create Differential Equation System

The standard way to create DE system is to create a derivative function in the form of f(du,u,p,t) with preallocated derivatives (du), initial variables ('u'), parameters (p) and time ('t')

function bistable_ode!(du, u, p ,t)
s1, s2 = u
K1, K2, k1, k2, k3, k4, n1 , n2  = p
du[1] = k1 / (1 + (s2/K2)^n1) - k3*s1
du[2] = k2/  (1 + (s1/K1)^n2) - k4*s2
end
NameTuple may cause error

One may want to use NameTuple like LabelledArrays to create Tuple array and apply to the ode function. However, this may cause getindex undefined error for Jacobian Calculation via ModelingToolkit, but still safe with solving steady states with NameTuple as returns.

## Create DE Problem for Solving

de = DEsteady(func = bistable_ode!,
p =  [1.,1.,20.,20.,5.,5.,4.,4.],
u0 = [3.,1.],
method = SSRootfind()
)

The DEsteady.method will further be used for DifferentialEquation.solve. Both dynamic solver (Dynamicss which reaches steady states by solving DE) and root finding (SSRootfind) can get the stable point. However, only SSRootfind can find the saddle point if the initial values is near by.

The system is now well organized and ready to be solved.

sol = solve(de)
u: 2-element Vector{Float64}:
3.999999765270408
0.015564205973794263

## Exploring the initial states

To find all the steady states, one needs to sample the initial state by grid or random search. To begin with, the ParameterGrid (grid search) and ParameterRandom (random search) are useful generator for iterating and returning the parameter set.

    param_rand = ParameterRandom(
methods= [
Uniform(0.,100.),
Log_uniform(0.001,100.)
],
len= 100
)
100-element ParameterRandom:
[24.35987673478841, 2.346518234169877]
[12.815682704126697, 1.9670007373318323]
[39.978536251948825, 99.07942522523452]
[31.482679067327844, 95.79460114896938]
[65.14362211440334, 0.018610679106675933]
[6.403263704023998, 0.004743911569750526]
[33.6421007178003, 2.0333940766390586]
[58.28925049633143, 1.114099465349585]
[24.35938669614075, 0.02999796439283893]
[8.37601852204144, 0.0012776428013927738]
⋮
[17.385926332034995, 0.0017145924864299248]
[32.098389458310095, 0.8833154996316388]
[0.30900110714018325, 0.021759874707872405]
[56.06934214848218, 0.11007644363961876]
[83.50494016483592, 72.93141530958678]
[73.80165362520121, 7.455134494126275]
[29.591636491270368, 0.26929879719681454]
[27.776464750956855, 2.8572021084758985]
[58.649895555775025, 0.0045185366363317166]
     param_grid = ParameterGrid([
[0.,100.,100],
[0.,100.,100]
])
10000-element ParameterGrid:
[0.0, 0.0]
[1.0101010101010102, 0.0]
[2.0202020202020203, 0.0]
[3.0303030303030303, 0.0]
[4.040404040404041, 0.0]
[5.050505050505051, 0.0]
[6.0606060606060606, 0.0]
[7.070707070707071, 0.0]
[8.080808080808081, 0.0]
[9.090909090909092, 0.0]
⋮
[91.91919191919193, 100.0]
[92.92929292929294, 100.0]
[93.93939393939395, 100.0]
[94.94949494949496, 100.0]
[95.95959595959597, 100.0]
[96.96969696969697, 100.0]
[97.97979797979798, 100.0]
[98.98989898989899, 100.0]
[100.0, 100.0]

    sols_rand = solve(de, param_rand)
EnsembleSolution Solution of length 100 with uType:
SciMLBase.NonlinearSolution{Float64, 1, Vector{Float64}, Vector{Float64}, SteadyStateProblem{Vector{Float64}, true, Vector{Float64}, ODEFunction{true, typeof(Main.bistable_ode!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, SSRootfind{SteadyStateDiffEq.var"#2#4"}, Nothing, Nothing}
    sols_grid = solve(de, param_grid)
EnsembleSolution Solution of length 10000 with uType:
SciMLBase.NonlinearSolution{Float64, 1, Vector{Float64}, Vector{Float64}, SteadyStateProblem{Vector{Float64}, true, Vector{Float64}, ODEFunction{true, typeof(Main.bistable_ode!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}}, SSRootfind{SteadyStateDiffEq.var"#2#4"}, Nothing, Nothing}