diff --git a/matlab/dyn_ramsey_static.m b/matlab/dyn_ramsey_static.m
index c2b619b3f..8f2ccea91 100644
--- a/matlab/dyn_ramsey_static.m
+++ b/matlab/dyn_ramsey_static.m
@@ -143,8 +143,8 @@ end
 % the auxiliary variables before the Lagrange multipliers are treated
 % as ordinary endogenous variables
 aux_eq = [1:orig_endo_aux_nbr, orig_endo_aux_nbr+orig_eq_nbr+1:size(fJ,1)];
-A = fJ(aux_eq,orig_endo_aux_nbr+1:end);
-y = res(aux_eq);
+A = fJ(1:orig_endo_nbr,orig_endo_nbr+find(aux_vars_type==6));
+y = res(1:orig_endo_nbr);
 mult = -A\y;
 
 resids1 = y+A*mult;
diff --git a/matlab/evaluate_steady_state.m b/matlab/evaluate_steady_state.m
index aec3a1b3f..6f240898b 100644
--- a/matlab/evaluate_steady_state.m
+++ b/matlab/evaluate_steady_state.m
@@ -125,7 +125,7 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
         %check whether steady state really solves the model
         resids = evaluate_static_model(ys,exo_ss,params,M,options);
 
-        n_multipliers=M.endo_nbr-M.orig_endo_nbr;
+        n_multipliers=M.orig_eq_nbr;
         nan_indices_multiplier=find(isnan(resids(1:n_multipliers)));
         nan_indices=find(isnan(resids(n_multipliers+1:end)));
 
@@ -196,9 +196,12 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
     elseif (options.bytecode == 0 && options.block == 0)
         if options.linear == 0
             % non linear model
-            [ys,check] = dynare_solve([M.fname '_static'],...
+            static_model = str2func([M.fname '_static']);
+            [ys,check] = dynare_solve(@static_problem,...
                                       ys_init,...
-                                      options, exo_ss, params);
+                                      options, exo_ss, params,...
+                                      M.orig_endo_nbr,...
+                                      static_model);
         else
             % linear model
             fh_static = str2func([M.fname '_static']);
@@ -294,3 +297,8 @@ function [ys,params,info] = evaluate_steady_state(ys_init,M,options,oo,steadysta
         info(2) = NaN;
         return
     end
+
+function [resids,jac] = static_problem(y,x,params,nvar,fh_static_model)
+    [r,j] = fh_static_model(y,x,params);
+    resids = r(1:nvar);
+    jac = j(1:nvar,1:nvar);
diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc
index 5a2fbdec7..bd14af7ce 100644
--- a/preprocessor/StaticModel.cc
+++ b/preprocessor/StaticModel.cc
@@ -1052,17 +1052,11 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
   initializeVariablesAndEquations();
 
   vector<BinaryOpNode *> neweqs;
-  for (unsigned int eq = 0; eq < equations.size(); eq++)
-    if (eq < equations.size() - aux_equations.size())
-      {
-        expr_t eq_tmp = equations[eq]->substituteStaticAuxiliaryVariable();
-        neweqs.push_back(dynamic_cast<BinaryOpNode *>(eq_tmp->toStatic(*this)));
-      }
-    else
-      {
-        expr_t eq_tmp = equations[eq]->substituteStaticAuxiliaryDefinition();
-        neweqs.push_back(dynamic_cast<BinaryOpNode *>(eq_tmp->toStatic(*this)));
-      }
+  for (unsigned int eq = 0; eq < equations.size() - aux_equations.size(); eq++)
+    {
+      expr_t eq_tmp = equations[eq]->substituteStaticAuxiliaryVariable();
+      neweqs.push_back(dynamic_cast<BinaryOpNode *>(eq_tmp->toStatic(*this)));
+    }
 
   equations.clear();
   copy(neweqs.begin(),neweqs.end(),back_inserter(equations));
@@ -1070,8 +1064,12 @@ StaticModel::computingPass(const eval_context_t &eval_context, bool no_tmp_terms
   set<int> vars;
 
   for (int i = 0; i < symbol_table.endo_nbr(); i++)
-    vars.insert(getDerivID(symbol_table.getID(eEndogenous, i), 0));
-
+    {
+      int id = symbol_table.getID(eEndogenous, i);
+      if (!symbol_table.isAuxiliaryVariableButNotMultiplier(id))
+                vars.insert(getDerivID(id, 0));
+    }        
+ 
   // Launch computations
   cout << "Computing static model derivatives:" << endl
        << " - order 1" << endl;
@@ -2126,7 +2124,7 @@ void StaticModel::writeAuxVarRecursiveDefinitions(const string &basename, const
 
   for (int i = 0; i < (int) aux_equations.size(); i++)
     {
-      dynamic_cast<ExprNode *>(aux_equations[i])->writeOutput(output, oMatlabStaticModel, temporary_terms, tef_terms);
+      dynamic_cast<ExprNode *>(aux_equations[i]->substituteStaticAuxiliaryDefinition())->writeOutput(output, oMatlabStaticModel, temporary_terms, tef_terms);
       output << ";" << endl;
     }
 }
diff --git a/preprocessor/SymbolTable.cc b/preprocessor/SymbolTable.cc
index 79fa54881..0a15fc070 100644
--- a/preprocessor/SymbolTable.cc
+++ b/preprocessor/SymbolTable.cc
@@ -786,6 +786,15 @@ SymbolTable::isAuxiliaryVariable(int symb_id) const
   return false;
 }
 
+bool
+SymbolTable::isAuxiliaryVariableButNotMultiplier(int symb_id) const
+{
+  for (int i = 0; i < (int) aux_vars.size(); i++)
+    if (aux_vars[i].get_symb_id() == symb_id && aux_vars[i].get_type() != avMultiplier)
+      return true;
+  return false;
+}
+
 set<int>
 SymbolTable::getOrigEndogenous() const
 {
diff --git a/preprocessor/SymbolTable.hh b/preprocessor/SymbolTable.hh
index 826d7a986..4f623dfdd 100644
--- a/preprocessor/SymbolTable.hh
+++ b/preprocessor/SymbolTable.hh
@@ -325,6 +325,8 @@ public:
   set <int> getEndogenous() const;
   //! Is a given symbol an auxiliary variable
   bool isAuxiliaryVariable(int symb_id) const;
+  //! Is a given symbol an auxiliary variable but not a Lagrange multiplier
+  bool isAuxiliaryVariableButNotMultiplier(int symb_id) const;
   //! Get list of endogenous variables without aux vars
   set <int> getOrigEndogenous() const;
 };