From 24d0f58d01b128cc40b41735ad999cc9e12c5b9e Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <Contact@ChrisRackauckas.com>
Date: Sat, 14 Apr 2018 10:55:53 -0700
Subject: [PATCH] add dynamicss

---
 src/SteadyStateDiffEq.jl |  6 +++---
 src/algorithms.jl        |  7 +++++++
 src/solve.jl             |  7 +++++++
 test/runtests.jl         | 12 ++++++++++++
 4 files changed, 29 insertions(+), 3 deletions(-)

diff --git a/src/SteadyStateDiffEq.jl b/src/SteadyStateDiffEq.jl
index 7316a69..1715a3c 100644
--- a/src/SteadyStateDiffEq.jl
+++ b/src/SteadyStateDiffEq.jl
@@ -4,8 +4,8 @@ module SteadyStateDiffEq
 
 using Reexport
 @reexport using DiffEqBase
-  
-using NLsolve
+
+using NLsolve, DiffEqCallbacks
 
 using Compat
 
@@ -14,6 +14,6 @@ import DiffEqBase: solve
 include("algorithms.jl")
 include("solve.jl")
 
-export SSRootfind
+export SSRootfind, DynamicSS
 
 end # module
diff --git a/src/algorithms.jl b/src/algorithms.jl
index 3fb057c..2832fd5 100644
--- a/src/algorithms.jl
+++ b/src/algorithms.jl
@@ -4,3 +4,10 @@ struct SSRootfind{F} <: SteadyStateDiffEqAlgorithm
   nlsolve::F
 end
 SSRootfind(;nlsolve=(f,u0,abstol) -> (res=NLsolve.nlsolve(f,u0,ftol = abstol);res.zero)) = SSRootfind(nlsolve)
+
+struct DynamicSS{A,AT,RT} <: SteadyStateDiffEqAlgorithm
+  alg::A
+  abstol::AT
+  reltol::RT
+end
+DynamicSS(alg;abstol = 1e-8, reltol = 1e-6) = DynamicSS(alg,abstol,reltol)
diff --git a/src/solve.jl b/src/solve.jl
index 6d36edb..64e645c 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -37,3 +37,10 @@ function solve(prob::AbstractSteadyStateProblem,alg::SteadyStateDiffEqAlgorithm,
     error("Algorithm not recognized")
   end
 end
+
+function solve(prob::AbstractSteadyStateProblem,alg::DynamicSS,args...;kwargs...)
+  _prob = ODEProblem(prob.f,prob.u0,(0.0,Inf),prob.p,
+                   mass_matrix=prob.mass_matrix,
+                   jac_prototype=prob.jac_prototype)
+  solve(_prob,alg.alg,args...;kwargs...,callback=TerminateSteadyState(alg.abstol,alg.reltol))
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index b8e99ca..f9be033 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -9,6 +9,7 @@ u0 = zeros(2)
 prob = SteadyStateProblem(f,u0)
 abstol = 1e-8
 sol = solve(prob,SSRootfind())
+p = nothing
 
 du = zeros(2)
 f(du,sol.u,p,0)
@@ -26,3 +27,14 @@ sol = solve(prob,SSRootfind(nlsolve = (f,u0,abstol) -> (res=Sundials.kinsol(f,u0
 
 f(du,sol.u,p,0)
 @test du == [0,0]
+
+using OrdinaryDiffEq
+sol = solve(prob,DynamicSS(Rodas5()))
+
+f(du,sol.u[end],p,0)
+@test du ≈ [0,0] atol = 1e-7
+
+sol = solve(prob,DynamicSS(CVODE_BDF()),dt=1.0)
+
+f(du,sol.u[end],p,0)
+@test du ≈ [0,0] atol = 1e-6