using ModelingToolkit
using SteadyStateDiffEq
using Test

@parameters t r
@variables x(t)
@derivatives D'~t
eqs = [D(x) ~ x^2-r]
de = ODESystem(eqs)

for factor in [1e-1, 1e0, 1e10], u0_p in [(2.34,2.676),(22.34,1.632),(.3,15.676),(0.3,0.006)]
    u0 = [x => factor*u0_p[1]]
    p = [r => factor*u0_p[2]]
    ss_prob = SteadyStateProblem(de,u0,p)
    sol = solve(ss_prob,SSRootfind()).u[1]
    @test abs(sol^2 - factor*u0_p[2]) < 1e-8
    ss_prob = SteadyStateProblemExpr(de,u0,p)
    sol_expr = solve(eval(ss_prob),SSRootfind()).u[1]
    @test all(x->x==0,sol-sol_expr)
end
