diff --git a/julia/DynareModel.jl b/julia/DynareModel.jl
index a85442c9d2bfe3741f7237d03720492e3b61077f..f89cb9effdfd3c55e38133d627ce7d859dc1a6d1 100644
--- a/julia/DynareModel.jl
+++ b/julia/DynareModel.jl
@@ -18,8 +18,7 @@ module DynareModel
  # along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 ##
 
-
-export Endo, Exo, ExoDet, Param, dynare_model
+export Model, Endo, Exo, ExoDet, Param, dynare_model
 
 abstract Atom
 
diff --git a/julia/DynareOptions.jl b/julia/DynareOptions.jl
index bdcca5c1621b3334526537040265705f88100d25..c5cc5e8881e91925f69b533e89ae2d1539b389c6 100644
--- a/julia/DynareOptions.jl
+++ b/julia/DynareOptions.jl
@@ -19,7 +19,7 @@ module DynareOptions
 ##
 
 
-export dynare_options
+export Options, dynare_options
 
 type Options
     dynare_version::ASCIIString
diff --git a/julia/DynareOutput.jl b/julia/DynareOutput.jl
index 81a4132f75017db245368b5fe00bc8904727f468..b3ca48a23a3ea399b792d37b881774376375f8e9 100644
--- a/julia/DynareOutput.jl
+++ b/julia/DynareOutput.jl
@@ -18,8 +18,7 @@ module DynareOutput
  # along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
 ##
 
-
-export dynare_output
+export Ouput, dynare_output
 
 type Output
     dynare_version::ASCIIString
@@ -28,7 +27,7 @@ type Output
 end
 
 function dynare_output()
-    return Output("",                   # dynare_version
+    return Output("",                # dynare_version
                   Array(Float64, 0), # steady_state
                   Array(Float64, 0)  # exo_steady_state
                  )
diff --git a/julia/SteadyState.jl b/julia/SteadyState.jl
new file mode 100644
index 0000000000000000000000000000000000000000..8d628563b6074ebbe5079189176e2bdd69ee567c
--- /dev/null
+++ b/julia/SteadyState.jl
@@ -0,0 +1,96 @@
+module SteadyState
+
+##
+ # Copyright (C) 2016 Dynare Team
+ #
+ # This file is part of Dynare.
+ #
+ # Dynare is free software: you can redistribute it and/or modify
+ # it under the terms of the GNU General Public License as published by
+ # the Free Software Foundation, either version 3 of the License, or
+ # (at your option) any later version.
+ #
+ # Dynare is distributed in the hope that it will be useful,
+ # but WITHOUT ANY WARRANTY; without even the implied warranty of
+ # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ # GNU General Public License for more details.
+ #
+ # You should have received a copy of the GNU General Public License
+ # along with Dynare.  If not, see <http://www.gnu.org/licenses/>.
+##
+
+using NLsolve
+
+import DynareModel.Model
+import DynareOutput.Output
+
+export steady, steady!
+export steady_state, steady_state!
+
+function steady(model::Model, oo::Output)
+    if model.analytical_steady_state || model.user_written_analytical_steady_state
+        steadystate = zeros(length(model.endo))
+        model.steady_state(steadystate, oo.exo_steady_state, model.params)
+        return steadystate
+    else
+        error("You have to provide a closed form solution for the steady state, or declare a guess\nfor the steady state as a third input argument.")
+    end
+end
+
+function steady!(model::Model, oo::Output)
+    if model.analytical_steady_state || model.user_written_analytical_steady_state
+        model.steady_state(oo.steady_state, oo.exo_steady_state, model.params)
+        return
+    else
+        error("You have to provide a closed form solution for the steady state, or declare a guess\nfor the steady state as a third input argument.")
+    end
+end
+
+function steady(model::Model, oo::Output, yinit::Vector{Float64})
+    ojectivefunction!(y::Vector{Float64}, fval::Vector{Float64}, fjac::Array{Float64}) = model.static(y, oo.exo_steady_state, model.params, fval, fjac)
+    r = nlsolve(only_fg!(ojectivefunction!), yinit, show_trace=false)
+    if converged(r)
+        return r.zero
+    else
+        return fill(nan(Float64), length(yinit))
+    end
+end
+
+function steady!(model::Model, oo::Output, yinit::Vector{Float64})
+    ojectivefunction!(y::Vector{Float64}, fval::Vector{Float64}, fjac::Array{Float64}) = model.static(y, oo.exo_steady_state, model.params, fval, fjac)
+    r = nlsolve(only_fg!(ojectivefunction!), yinit, show_trace=false)
+    if converged(r)
+        oo.steady_state = r.zero
+    else
+        oo.steady_state = fill(nan(Float64), length(yinit))
+    end
+end
+
+function steady_state(model::Model, oo::Output)
+    ys = steady(model, oo)
+    display_steady_state(model, oo, ys)
+end
+
+function steady_state!(model::Model, oo::Output)
+    steady!(model, oo)
+    display_steady_state(model, oo, oo.steady_state)
+end
+
+function display_steady_state(model::Model, oo::Output, ys::Vector{Float64})
+    println("\n\nSTEADY STATE:\n")
+    for i = 1:length(model.endo)
+        println(string(model.endo[i].name,  " = ",  ys[i]))
+    end
+end
+
+function issteadystate(model::Model, oo::Output, ys::Vector{Float64})
+    residuals = zeros(Float64, length(ys))
+    compute_static_model_residuals!(model, oo, ys, residuals)
+    return maximum(abs(residuals))<1e-6
+end
+
+function compute_static_model_residuals!(model::Model, oo::Output, ys::Vector{Float64}, residuals::Vector{Float64})
+    model.static(ys, oo.exo_steady_state, model.params, residuals)
+end
+
+end
diff --git a/preprocessor/ModFile.cc b/preprocessor/ModFile.cc
index 2b2c684bdb2ef67b50541ff81aa487bb135a3216..05c7bce25fd712ebb938de01abe51ba2b92fc496 100644
--- a/preprocessor/ModFile.cc
+++ b/preprocessor/ModFile.cc
@@ -1104,11 +1104,12 @@ ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output)
                << "#" << endl
                << "# NB: this file was automatically generated by Dynare" << endl
                << "#     from " << basename << ".mod" << endl
-               << "#" << endl
+               << "#" << endl << endl
                << "using DynareModel" << endl
                << "using DynareOptions" << endl
-               << "using DynareOutput" << endl
+               << "using DynareOutput" << endl << endl
                << "using Utils" << endl
+               << "using SteadyState" << endl << endl
                << "using " << basename << "Static" << endl
                << "using " << basename << "Dynamic" << endl
                << "if isfile(\"" << basename << "SteadyState.jl"  "\")" << endl
@@ -1197,7 +1198,7 @@ ModFile::writeExternalFilesJulia(const string &basename, FileOutputType output)
                << "    using " << basename << "DynamicParamsDerivs" << endl
                << "    model_.dynamic_params_derivs = " << basename << "DynamicParamsDerivs.params_derivs" << endl
                << "end" << endl
-               << "end" << endl;
+	       << "end" << endl;
   jlOutputFile.close();
   cout << "done" << endl;
 }
diff --git a/tests/julia/rbc/test1.jl b/tests/julia/rbc/test1.jl
index 8f7ba07ca560b2fbb95769b6a9a016e3909103b4..c4579fa377de887f1603f09081abde5ba406ee3f 100644
--- a/tests/julia/rbc/test1.jl
+++ b/tests/julia/rbc/test1.jl
@@ -10,4 +10,83 @@ importall Dynare
 
 @dynare "rbc1.mod"
 
-print(model_.fname)
+importall SteadyState
+
+# First call to the steady state routine (analytical)
+@time SteadyState.steady!(model_, oo_)
+
+# First call to the steady state routine (analytical)
+@time SteadyState.steady!(model_, oo_)
+
+paramsinit = copy(model_.params);
+
+yinit = deepcopy(oo_.steady_state)
+
+y_init = copy(yinit)
+y_init[1] = 1.1*yinit[1]
+y_init[2] = 0.9*yinit[2]
+
+# First call to the steady state routine (numerical)
+println("First call to the numerical steady state routine")
+@time SteadyState.steady!(model_, oo_, yinit)
+
+# Check results
+@assert maximum(abs(oo_.steady_state-yinit))<1e-9
+
+y_init = copy(yinit)
+yinit[1] = 1.1*yinit[1]
+yinit[2] = 0.9*yinit[2]
+
+# Second call to the steady state routine (numerical)
+println("Second call to the numerical steady state routine")
+@time SteadyState.steady!(model_, oo_, yinit)
+
+# change alpha
+println("Change α")
+model_.params[4] = max(min(1.0, model_.params[4]*1.2), 0.0)
+ys = SteadyState.steady(model_, oo_)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+@assert maximum(abs(oo_.steady_state-ys))<1e-9
+
+# change delta
+println("Change δ")
+model_.params[6] = max(min(1.0, model_.params[6]*1.2), 0.0)
+ys = SteadyState.steady(model_, oo_)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+@assert maximum(abs(oo_.steady_state-ys))<1e-9
+
+# change beta
+println("Change β")
+model_.params[1] = max(min(1-1e-6, model_.params[1]*0.99), 0.0)
+ys = SteadyState.steady(model_, oo_)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+@assert maximum(abs(oo_.steady_state-ys))<1e-9
+
+# change tau
+println("Change τ")
+model_.params[3] /= 1.5
+ys = SteadyState.steady(model_, oo_)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+@assert maximum(abs(oo_.steady_state-ys))<1e-9
+
+# change Epsilon
+println("Change ϵ")
+model_.params[5] *= 1.5
+ys = SteadyState.steady(model_, oo_)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+y_init = copy(yinit)
+@time SteadyState.steady!(model_, oo_, y_init)
+@assert maximum(abs(oo_.steady_state-ys))<1e-9
diff --git a/tests/julia/rbc/test2.jl b/tests/julia/rbc/test2.jl
index aa7a9566bd952e432258cc9420daf4a985c32afc..f9eae9ebfad5e472ca4890461f6395484a7503cd 100644
--- a/tests/julia/rbc/test2.jl
+++ b/tests/julia/rbc/test2.jl
@@ -11,4 +11,6 @@ importall Dynare
 
 @dynare "rbc2.mod"
 
-print(model_.fname)
+importall SteadyState
+
+steady!(model_, oo_)