diff --git a/matlab/dynare_m.exe b/matlab/dynare_m.exe
index acb97b6ddc3d14bcfb9719cc7cee9b2bdf83aa0d..101ca309ad7d5b1ac3ffecf4978e0834e4639444 100755
Binary files a/matlab/dynare_m.exe and b/matlab/dynare_m.exe differ
diff --git a/matlab/homotopy1.m b/matlab/homotopy1.m
index 3e928968c8f6f4c8c8f1f1b45011edd422b04595..fb38691b6b63165fe1c8f9fc191790da2b51cb48 100644
--- a/matlab/homotopy1.m
+++ b/matlab/homotopy1.m
@@ -13,6 +13,8 @@ function homotopy1(values, step_nbr)
 %                   exogenous deterministic, 4 for parameters)
 %                   Column 2 is symbol integer identifier.
 %                   Column 3 is initial value, and column 4 is final value.
+%                   Column 3 can contain NaNs, in which case previous
+%                   initialization of variable will be used as initial value.
 %    step_nbr:      number of steps for homotopy
 %
 % OUTPUTS
@@ -36,15 +38,25 @@ function homotopy1(values, step_nbr)
     error('HOMOTOPY: incorrect variable types specified')
   end
 
-  if any(values(:,3) == values(:,4))
+  % Construct vector of starting values, using previously initialized values
+  % when initial value has not been given in homotopy_setup block
+  oldvalues = values(:,3);
+  ipn = find(values(:,1) == 4 & isnan(oldvalues));
+  oldvalues(ipn) = M_.params(values(ipn, 2));
+  ixn = find(values(:,1) == 1 & isnan(oldvalues));
+  oldvalues(ixn) = oo_.exo_steady_state(values(ixn, 2));
+  ixdn = find(values(:,1) == 2 & isnan(oldvalues));
+  oldvalues(ixdn) = oo_.exo_det_steady_state(values(ixdn, 2));
+  
+  if any(oldvalues == values(:,4))
     error('HOMOTOPY: initial and final values should be different')
   end
   
   points = zeros(nv, step_nbr+1);
   for i = 1:nv
-    points(i,:) = values(i,3):(values(i,4)-values(i,3))/step_nbr:values(i,4);
+    points(i,:) = oldvalues(i):(values(i,4)-oldvalues(i))/step_nbr:values(i,4);
   end
-    
+  
   for i=1:step_nbr+1
     M_.params(values(ip,2)) = points(ip,i);
     oo_.exo_steady_state(values(ix,2)) = points(ix,i);
diff --git a/matlab/homotopy2.m b/matlab/homotopy2.m
index 3518a2a71fe2e416fcc7111075f77fc9d06e88d4..813d8c1465370f45b082604dc5b21d0475772672 100755
--- a/matlab/homotopy2.m
+++ b/matlab/homotopy2.m
@@ -15,6 +15,8 @@ function homotopy2(values, step_nbr)
 %                   exogenous deterministic, 4 for parameters)
 %                   Column 2 is symbol integer identifier.
 %                   Column 3 is initial value, and column 4 is final value.
+%                   Column 3 can contain NaNs, in which case previous
+%                   initialization of variable will be used as initial value.
 %    step_nbr:      number of steps for homotopy
 %
 % OUTPUTS
@@ -30,27 +32,41 @@ function homotopy2(values, step_nbr)
 
   nv = size(values, 1);
   
-  % Initialize all variables with initial value
+  oldvalues = values(:,3);
+  
+  % Initialize all variables with initial value, or the reverse...
   for i = 1:nv
     switch values(i,1)
      case 1
-      oo_.exo_steady_state(values(i,2)) = values(i,3);
+      if isnan(oldvalues(i))
+        oldvalues(i) = oo_.exo_steady_state(values(i,2));
+      else
+        oo_.exo_steady_state(values(i,2)) = oldvalues(i);
+      end
      case 2
-      oo_.exo_det_steady_state(values(i,2)) = values(i,3);
+      if isnan(oldvalues(i))
+        oldvalues(i) = oo_.exo_det_steady_state(values(i,2));
+      else
+        oo_.exo_det_steady_state(values(i,2)) = oldvalues(i);
+      end
      case 4
-      M_.params(values(i,2)) = values(i,3);
+      if isnan(oldvalues(i))
+        oldvalues(i) = M_.params(values(i,2));
+      else
+        M_.params(values(i,2)) = oldvalues(i);
+      end
      otherwise
       error('HOMOTOPY: incorrect variable types specified')
     end
   end
 
-  if any(values(:,3) == values(:,4))
+  if any(oldvalues == values(:,4))
     error('HOMOTOPY: initial and final values should be different')
   end
   
   % Actually do the homotopy
   for i = 1:nv
-    for v = values(i,3):(values(i,4)-values(i,3))/step_nbr:values(i,4)
+    for v = oldvalues(i):(values(i,4)-oldvalues(i))/step_nbr:values(i,4)
       switch values(i,1)
        case 1
         oo_.exo_steady_state(values(i,2)) = v;
diff --git a/matlab/homotopy3.m b/matlab/homotopy3.m
index d53c1b1b99be0fc962fbdebc77e6569799781bd6..f895970ca24d2faf688e36ca171da596d1b1a871 100644
--- a/matlab/homotopy3.m
+++ b/matlab/homotopy3.m
@@ -15,6 +15,8 @@ function homotopy3(values, step_nbr)
 %                   exogenous deterministic, 4 for parameters)
 %                   Column 2 is symbol integer identifier.
 %                   Column 3 is initial value, and column 4 is final value.
+%                   Column 3 can contain NaNs, in which case previous
+%                   initialization of variable will be used as initial value.
 %    step_nbr:      maximum number of steps to try before aborting
 %
 % OUTPUTS
@@ -39,7 +41,16 @@ function homotopy3(values, step_nbr)
     error('HOMOTOPY: incorrect variable types specified')
   end
 
+  % Construct vector of starting values, using previously initialized values
+  % when initial value has not been given in homotopy_setup block
   oldvalues = values(:,3);
+  ipn = find(values(:,1) == 4 & isnan(oldvalues));
+  oldvalues(ipn) = M_.params(values(ipn, 2));
+  ixn = find(values(:,1) == 1 & isnan(oldvalues));
+  oldvalues(ixn) = oo_.exo_steady_state(values(ixn, 2));
+  ixdn = find(values(:,1) == 2 & isnan(oldvalues));
+  oldvalues(ixdn) = oo_.exo_det_steady_state(values(ixdn, 2));
+
   targetvalues = values(:,4);
 
   if min(abs(targetvalues - oldvalues)) < tol
diff --git a/preprocessor/DynareBison.yy b/preprocessor/DynareBison.yy
index 8ee0389b57ac5cfa7723f92385748e3fe16c42d0..df1d842400b9c5d5cbd32ad6c0e7789d03cd14a7 100644
--- a/preprocessor/DynareBison.yy
+++ b/preprocessor/DynareBison.yy
@@ -1250,7 +1250,10 @@ homotopy_list : homotopy_item
               ;
 
 homotopy_item : NAME COMMA expression COMMA expression ';'
-              { driver.homotopy_val($1,$3,$5);};
+                { driver.homotopy_val($1, $3, $5);}
+              | NAME COMMA expression ';'
+                { driver.homotopy_val($1, NULL, $3);}
+              ;
 
 number : INT_NUMBER
        | FLOAT_NUMBER
diff --git a/preprocessor/NumericalInitialization.cc b/preprocessor/NumericalInitialization.cc
index b7783141e8130aff7eb084b091d5c9e3967f3cbf..a38f4fdaffb56064b4e2a7ed3ddf20e7bd04d7d9 100644
--- a/preprocessor/NumericalInitialization.cc
+++ b/preprocessor/NumericalInitialization.cc
@@ -179,7 +179,10 @@ HomotopyStatement::writeOutput(ostream &output, const string &basename) const
       const int id = symbol_table.getID(name) + 1;
 
       output << "options_.homotopy_values = vertcat(options_.homotopy_values, [ " << type << ", " << id << ", ";
-      expression1->writeOutput(output);
+      if (expression1 != NULL)
+        expression1->writeOutput(output);
+      else
+        output << "NaN";
       output << ", ";
       expression2->writeOutput(output);
       output << "]);" << endl;
diff --git a/preprocessor/ParsingDriver.cc b/preprocessor/ParsingDriver.cc
index eec83386b03ce380ba4e9d160588cb1db1f0382d..6e0e10b120ad487d4bd8f3cd8a8d0c71082f58f9 100644
--- a/preprocessor/ParsingDriver.cc
+++ b/preprocessor/ParsingDriver.cc
@@ -347,10 +347,7 @@ ParsingDriver::homotopy_val(string *name, NodeID val1, NodeID val2)
       && type != eExogenousDet)
     error("homotopy_val: " + *name + " should be a parameter or exogenous variable");
 
-  if (homotopy_values.find(*name) != homotopy_values.end())
-    error("homotopy_val: " + *name +" declared twice");
-
-  homotopy_values[*name] = make_pair(val1, val2);
+  homotopy_values.push_back(make_pair(*name, make_pair(val1, val2)));
 
   delete name;
 }
diff --git a/preprocessor/include/NumericalInitialization.hh b/preprocessor/include/NumericalInitialization.hh
index f920c96756df6553ec6da7a70808361b102f65be..41459423a8d399e880cb9b4029d57b0e82283231 100644
--- a/preprocessor/include/NumericalInitialization.hh
+++ b/preprocessor/include/NumericalInitialization.hh
@@ -95,11 +95,9 @@ public:
 class HomotopyStatement : public Statement
 {
 public:
-  /*!
-    Contrary to Initval and Endval, we use a map, since it is impossible to reuse
-    a given initialization value in a second initialization inside the block.
-  */
-  typedef map<string,pair<NodeID,NodeID> > homotopy_values_type;
+  //! Stores the declarations of homotopy_setup
+  /*! Order matter so we use a vector. First NodeID can be NULL if no initial value given. */
+  typedef vector<pair<string, pair<NodeID, NodeID> > > homotopy_values_type;
 private:
   const homotopy_values_type homotopy_values;
   const SymbolTable &symbol_table;
diff --git a/preprocessor/include/ParsingDriver.hh b/preprocessor/include/ParsingDriver.hh
index 72c0cdbea71e045c6f8fb12afbc8a827ef9d1af5..071e5592eb5d0f91324b1a695827c24fd05878b2 100644
--- a/preprocessor/include/ParsingDriver.hh
+++ b/preprocessor/include/ParsingDriver.hh
@@ -209,7 +209,8 @@ public:
   void init_val(string *name, NodeID rhs);
   //! Writes an histval block
   void hist_val(string *name, string *lag, NodeID rhs);
-  //! Writes an homotopy_setup block
+  //! Adds an entry in a homotopy_setup block
+  /*! Second argument "val1" can be NULL if no initial value provided */
   void homotopy_val(string *name, NodeID val1, NodeID val2);
   //! Writes end of an initval block
   void end_initval();
diff --git a/tests/homotopy/ramst_homotopy.mod b/tests/homotopy/ramst_homotopy.mod
index 55bac89137eb1b086b7bf71447e35a4a694e137a..8b468f4495efc0755ff4095ec8de7f028026fcec 100755
--- a/tests/homotopy/ramst_homotopy.mod
+++ b/tests/homotopy/ramst_homotopy.mod
@@ -7,13 +7,14 @@
 initval;
 k = 12.75;
 c = 1.5;
+x = 1;
 end;
 
 homotopy_setup;
 bet, 0.05, 0.1;
-x, 1, 2;
+x, 2;
 end;
 
-steady(homotopy_mode = 1, homotopy_steps = 50);
+steady(homotopy_mode = 2, homotopy_steps = 50);
 //steady(homotopy_mode = 2, homotopy_steps = 50);
 //steady(homotopy_mode = 3, homotopy_steps = 50);