diff --git a/doc/manual.xml b/doc/manual.xml
index 06f426a3a2a07f913e15b1118a0a2de72a2f05c5..d19391db84238faff3938e22d347574159fa9d98 100644
--- a/doc/manual.xml
+++ b/doc/manual.xml
@@ -1836,7 +1836,7 @@ steady(homotopy_mode = 1, homotopy_steps = 50);
       <itemizedlist>
         <listitem><para><literal>0</literal>: Newton method to solve simultaneously all the equations for every period, see <xref linkend="juillard_1996"/>. (Default)</para></listitem>
         <listitem><para><literal>1</literal>: use a Newton algorithm with a sparse LU solver at each iteration.</para></listitem>
-        <listitem><para><literal>2</literal>: use a Newton algorithm with a Generalized Minimal Residual (GMRES) solver at each iteration.</para></listitem>
+        <listitem><para><literal>2</literal>: use a Newton algorithm with a Generalized Minimal Residual (GMRES) solver at each iteration. This option is not available under Octave.</para></listitem>
         <listitem><para><literal>3</literal>: use a Newton algorithm with a Stabilized Bi-Conjugate Gradient (BICGSTAB) solver at each iteration.</para></listitem>
         <listitem><para><literal>4</literal>: use a Newton algorithm with a optimal path length at each iteration.</para></listitem>
         <listitem><para><literal>5</literal>: use a Newton algorithm with a sparse Gaussian elimination (SPE) solver at each iteration.</para></listitem>
diff --git a/matlab/lnsrch1_wrapper_one_boundary.m b/matlab/lnsrch1_wrapper_one_boundary.m
new file mode 100644
index 0000000000000000000000000000000000000000..fcc5c465a09fc4be45ad9d8586bf23327dbcb164
--- /dev/null
+++ b/matlab/lnsrch1_wrapper_one_boundary.m
@@ -0,0 +1,43 @@
+function r = lnsrch1_wrapper_one_boundary(ya, y_index, fname, y, x, params, it_)
+% wrapper for solve_one_boundary m-file when it is used with a dynamic
+% model
+%
+% INPUTS
+%   ya                  [vector]        The endogenous of the current block
+%   y_index             [vector of int] The index of the endogenous variables of
+%                                       the block
+%   fname               [string]        name of the file containing the block
+%                                       to simulate
+%   y                   [matrix]        All the endogenous variables of the model
+%   x                   [matrix]        All the exogenous variables of the model
+%   params              [vector]        All the parameters of the model
+% OUTPUTS
+%   r                   [vector]        The residuals of the current block      
+%  
+% ALGORITHM
+%   none.
+%    
+% SPECIAL REQUIREMENTS
+%   none.
+%  
+
+% Copyright (C) 2009 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/licen
+
+%reshape the input arguments of the dynamic function
+y(it_, :) = ya;
+[r, y, g1, g2, g3]=feval(fname, y, x, params, it_, 0);
diff --git a/matlab/lnsrch1_wrapper_two_boundaries.m b/matlab/lnsrch1_wrapper_two_boundaries.m
new file mode 100644
index 0000000000000000000000000000000000000000..d8e946cff645f7391f5c96e7a93672e6fba4cb36
--- /dev/null
+++ b/matlab/lnsrch1_wrapper_two_boundaries.m
@@ -0,0 +1,48 @@
+function ra = lnsrch1_wrapper_two_boundaries(ya, fname, y, y_index, x, params, periods, y_kmin, y_size)
+% wrapper for solve_one_boundary m-file when it is used with a dynamic
+% model
+%
+% INPUTS
+%   ya                  [vector]        The endogenous of the current block
+%   y_index             [vector of int] The index of the endogenous variables of
+%                                       the block
+%   fname               [string]        name of the file containing the block
+%                                       to simulate
+%   y                   [matrix]        All the endogenous variables of the model
+%   x                   [matrix]        All the exogenous variables of the model
+%   params              [vector]        All the parameters of the model
+%   periods             [int]           The number of periods
+%   y_kmin              [int]           The maximum number of lag on en endogenous variables
+%   y_size              [int]           The number of endogenous variables
+%                                       in the current block
+% OUTPUTS
+%   ra                  [vector]        The residuals of the current block      
+%  
+% ALGORITHM
+%   none.
+%    
+% SPECIAL REQUIREMENTS
+%   none.
+%  
+
+% Copyright (C) 2009 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/licen
+
+%reshape the input arguments of the dynamic function
+y(y_kmin+1:y_kmin+periods, y_index) = reshape(ya',length(y_index),periods)';
+[r, y, g1, g2, g3, b]=feval(fname, y, x, params, periods, 0, y_kmin, y_size);
+ra = reshape(r(:, y_kmin+1:periods+y_kmin),periods*y_size, 1);
diff --git a/matlab/simul.m b/matlab/simul.m
index b2d755881dcb5c55baa3fc22ca788fad661ed347..0f3e4f5312421346741e2ecffa8305e092e513a1 100644
--- a/matlab/simul.m
+++ b/matlab/simul.m
@@ -67,6 +67,10 @@ if options_.block && options_.bytecode && options_.stack_solve_algo ~= 5
     error('SIMUL: for the moment, you must use stack_solve_algo=5 with block and bytecode option')
 end
 
+if exist('OCTAVE_VERSION') && options_.stack_solve_algo == 2
+    error('SIMUL: stack_solve_algo=2 is not available for Octave. Choose another value.')
+end
+
 if(options_.block)
     if(options_.bytecode)
         oo_.endo_simul=bytecode('dynamic');
diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m
index d694f095829933e1b691dec26143129986b0cd26..90f5090a09c5a00439a55d14c4fbeed16c7f06e5 100644
--- a/matlab/solve_one_boundary.m
+++ b/matlab/solve_one_boundary.m
@@ -237,20 +237,25 @@ for it_=start:incr:finish
                 lambda=1;
                 stpmx = 100 ;
                 if (is_dynamic)
-                    stpmax = stpmx*max([sqrt(y*y');size(y_index_eq,2)]);
+                    stpmax = stpmx*max([sqrt(ya'*ya);size(y_index_eq,2)]);
                 else
-                    stpmax = stpmx*max([sqrt(y'*y);size(y_index_eq,2)]);
+                    stpmax = stpmx*max([sqrt(ya'*ya);size(y_index_eq,2)]);
                 end;
                 nn=1:size(y_index_eq,2);
                 g = (r'*g1)';
                 f = 0.5*r'*r;
                 p = -g1\r ;
                 if (is_dynamic)
-                    [y,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq,x, params, it_, 0);
+                    [ya,f,r,check]=lnsrch1(y(it_,:),f,g,p,stpmax,'lnsrch1_wrapper_one_boundary',nn,  y_index_eq, y_index_eq, fname, y, x, params, it_);
                 else
-                    [y,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq,x, params, 0);
+                    [ya,f,r,check]=lnsrch1(y,f,g,p,stpmax,fname,nn,y_index_eq,x, params, 0);
                 end;
                 dx = ya - y(y_index_eq);
+                if(is_dynamic)
+                    y(it_,:) = ya';
+                else
+                    y = ya';
+                end;
             elseif(~is_dynamic & options_.solve_algo==3)
                 [yn,info] = csolve(@local_fname, y(y_index_eq),@local_fname,1e-6,500, x, params, y, y_index_eq, fname, 1);
                 dx = ya - yn;
@@ -267,8 +272,7 @@ for it_=start:incr:finish
                 flag1=1;
                 while(flag1>0)
                     [L1, U1]=luinc(g1,luinc_tol);
-                    [za,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1);
-                    %[za,flag1] = gmres(-g1,b',Blck_size,1e-6,Blck_size,L1,U1);
+                    [dx,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1);
                     if (flag1>0 | reduced)
                         if(flag1==1)
                             disp(['Error in simul: No convergence inside GMRES after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')]);
@@ -280,7 +284,6 @@ for it_=start:incr:finish
                         luinc_tol = luinc_tol/10;
                         reduced = 0;
                     else
-                        dx = za - ya;
                         ya = ya + lambda*dx;
                         if(is_dynamic)
                             y(it_,y_index_eq) = ya';
@@ -293,8 +296,7 @@ for it_=start:incr:finish
                 flag1=1;
                 while(flag1>0)
                     [L1, U1]=luinc(g1,luinc_tol);
-                    [za,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1);
-                    %[za,flag1] = bicgstab(-g1,b',1e-7,Blck_size,L1,U1);
+                    [dx,flag1] = bicgstab(g1,-r,1e-7,Blck_size,L1,U1);
                     if (flag1>0 | reduced)
                         if(flag1==1)
                             disp(['Error in simul: No convergence inside BICGSTAB after ' num2str(iter,'%6d') ' iterations, in block' num2str(Block_Num,'%3d')]);
@@ -306,7 +308,6 @@ for it_=start:incr:finish
                         luinc_tol = luinc_tol/10;
                         reduced = 0;
                     else
-                        dx = za - ya;
                         ya = ya + lambda*dx;
                         if(is_dynamic)
                             y(it_,y_index_eq) = ya';
diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m
index 9b54d530cf46db44f9cdfd8b00eb1204c29aec9b..fa7dc79a2bfcc5affce5a098ab8f67119f30aac8 100644
--- a/matlab/solve_two_boundaries.m
+++ b/matlab/solve_two_boundaries.m
@@ -213,14 +213,6 @@ while ~(cvg==1 | iter>maxit_),
             dx = g1a\b- ya;
             ya = ya + lambda*dx;
             y(1+y_kmin:periods+y_kmin,y_index)=reshape(ya',length(y_index),periods)';
-            %         v = '';
-            %         for i=1:(size(y_index,2))
-            %           v = [v ' %1.6f'];
-            %         end;
-            %         v = [v '\n'];
-            %         v
-            %         sprintf(v,y(:,y_index)')
-            %         return;
         elseif(stack_solve_algo==2),
             flag1=1;
             while(flag1>0)
@@ -264,11 +256,20 @@ while ~(cvg==1 | iter>maxit_),
                 end;
             end;
         elseif(stack_solve_algo==4),
-            error('SOLVE_TWO_BOUNDARIES: stack_solve_algo=4 not implemented')
-        end;
+            ra = reshape(r(:, y_kmin+1:periods+y_kmin),periods*Blck_size, 1);
+            stpmx = 100 ;
+            stpmax = stpmx*max([sqrt(ya'*ya);size(y_index,2)]);
+            nn=1:size(ra,1);
+            g = (ra'*g1a)';
+            f = 0.5*ra'*ra;
+            p = -g1a\ra;
+            [yn,f,ra,check]=lnsrch1(ya,f,g,p,stpmax,'lnsrch1_wrapper_two_boundaries',nn,nn,  fname, y, y_index, x, params, periods, y_kmin, Blck_size);
+            dx = ya - yn;
+            y(1+y_kmin:periods+y_kmin,y_index)=reshape(yn',length(y_index),periods)';
+        end
     end
     iter=iter+1;
-    disp(['iteration: ' num2str(iter,'%d') ' error: ' num2str(max_res,'%e')]);
+    disp(['iteration: ' num2str(iter,'%d') ' error: ' num2str(max_res,'%e') ' stack_solve_algo=' num2str(stack_solve_algo)]);
 end;
 if (iter>maxit_)
     disp(['No convergence after ' num2str(iter,'%4d') ' iterations in Block ' num2str(Block_Num,'%d')]);
diff --git a/preprocessor/DynamicModel.cc b/preprocessor/DynamicModel.cc
index 26b7c2c91a08d5de8ea87755504f5620a4bbc308..1bfb09d63a578b03c7ff4f90d2c6077e64c1ac2e 100644
--- a/preprocessor/DynamicModel.cc
+++ b/preprocessor/DynamicModel.cc
@@ -1995,7 +1995,7 @@ DynamicModel::computingPass(bool jacobianExo, bool hessian, bool thirdDerivative
 
       evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false);
 
-      computePossiblySingularNormalization(contemporaneous_jacobian, cutoff == 0);
+      computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian);
 
       computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered, prologue, epilogue);
 
diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc
index 0b2853159de51a9d8820df22d081d51b413cc1bd..d7d6f1e28f786fe16566a561d8af992ddd77c0f0 100644
--- a/preprocessor/ModelTree.cc
+++ b/preprocessor/ModelTree.cc
@@ -32,8 +32,8 @@
 using namespace boost;
 using namespace MFS;
 
-void
-ModelTree::computeNormalization(const set<pair<int, int> > &endo_eqs_incidence) throw (NormalizationException)
+bool
+ModelTree::computeNormalization(const jacob_map &contemporaneous_jacobian, bool verbose)
 {
   const int n = equation_number();
 
@@ -50,8 +50,8 @@ ModelTree::computeNormalization(const set<pair<int, int> > &endo_eqs_incidence)
   // Fill in the graph
   set<pair<int, int> > endo;
 
-  for (set<pair<int, int> >::const_iterator it = endo_eqs_incidence.begin(); it != endo_eqs_incidence.end(); it++)
-    add_edge(it->first + n, it->second, g);
+  for (jacob_map::const_iterator it = contemporaneous_jacobian.begin(); it != contemporaneous_jacobian.end(); it++)
+    add_edge(it->first.first + n, it->first.second, g);
 
   // Compute maximum cardinality matching
   vector<int> mate_map(2*n);
@@ -125,63 +125,96 @@ ModelTree::computeNormalization(const set<pair<int, int> > &endo_eqs_incidence)
   // Check if all variables are normalized
   vector<int>::const_iterator it = find(mate_map.begin(), mate_map.begin() + n, graph_traits<BipartiteGraph>::null_vertex());
   if (it != mate_map.begin() + n)
-    throw NormalizationException(symbol_table.getID(eEndogenous, it - mate_map.begin()));
+    {
+      if (verbose)
+        cerr << "ERROR: Could not normalize the model. Variable "
+             << symbol_table.getName(symbol_table.getID(eEndogenous, it - mate_map.begin()))
+             << " is not in the maximum cardinality matching." << endl;
+      check = false;
+    }
+  return check;
 }
 
 void
-ModelTree::computePossiblySingularNormalization(const jacob_map &contemporaneous_jacobian, bool try_symbolic)
+ModelTree::computeNonSingularNormalization(jacob_map &contemporaneous_jacobian, double cutoff, jacob_map &static_jacobian, dynamic_jacob_map &dynamic_jacobian)
 {
+  bool check = false;
+
   cout << "Normalizing the model..." << endl;
 
-  set<pair<int, int> > endo_eqs_incidence;
+  int n = equation_number();
 
-  for (jacob_map::const_iterator it = contemporaneous_jacobian.begin();
-       it != contemporaneous_jacobian.end(); it++)
-    endo_eqs_incidence.insert(make_pair(it->first.first, it->first.second));
+  // compute the maximum value of each row of the contemporaneous Jacobian matrix
+  //jacob_map normalized_contemporaneous_jacobian;
+  jacob_map normalized_contemporaneous_jacobian(contemporaneous_jacobian);
+  vector<double> max_val(n, 0.0);
+  for (jacob_map::const_iterator iter = contemporaneous_jacobian.begin(); iter != contemporaneous_jacobian.end(); iter++)
+    if (fabs(iter->second) > max_val[iter->first.first])
+      max_val[iter->first.first] = fabs(iter->second);
 
-  try
-    {
-      computeNormalization(endo_eqs_incidence);
-      return;
-    }
-  catch (NormalizationException &e)
+  for (jacob_map::iterator iter = normalized_contemporaneous_jacobian.begin(); iter != normalized_contemporaneous_jacobian.end(); iter++)
+    iter->second /= max_val[iter->first.first];
+
+  //We start with the highest value of the cutoff and try to normalize the model
+  double current_cutoff = 0.99999999;
+
+  int suppressed = 0;
+  while (!check && current_cutoff > 1e-19)
     {
-      if (try_symbolic)
-        cout << "Normalization failed with cutoff, trying symbolic normalization..." << endl;
-      else
+      jacob_map tmp_normalized_contemporaneous_jacobian;
+      int suppress = 0;
+      for (jacob_map::iterator iter = normalized_contemporaneous_jacobian.begin(); iter != normalized_contemporaneous_jacobian.end(); iter++)
+        if (fabs(iter->second) > max(current_cutoff, cutoff))
+          tmp_normalized_contemporaneous_jacobian[make_pair(iter->first.first, iter->first.second)] = iter->second;
+        else
+          suppress++;
+
+      if (suppress != suppressed)
+        check = computeNormalization(tmp_normalized_contemporaneous_jacobian, false);
+      suppressed = suppress;
+      if (!check)
         {
-          cerr << "ERROR: Could not normalize the model. Variable "
-               << symbol_table.getName(e.symb_id)
-               << " is not in the maximum cardinality matching. Try to decrease the cutoff." << endl;
-          exit(EXIT_FAILURE);
+          current_cutoff /= 2;
+          // In this last case try to normalize with the complete jacobian
+          if (current_cutoff <= 1e-19)
+            check = computeNormalization(normalized_contemporaneous_jacobian, false);
         }
     }
 
-  // If no non-singular normalization can be found, try to find a normalization even with a potential singularity
-  if (try_symbolic)
+  if (!check)
     {
-      endo_eqs_incidence.clear();
+      cout << "Normalization failed with cutoff, trying symbolic normalization..." << endl;
+      //if no non-singular normalization can be found, try to find a normalization even with a potential singularity
+      jacob_map tmp_normalized_contemporaneous_jacobian;
       set<pair<int, int> > endo;
-      for (int i = 0; i < equation_number(); i++)
+      for (int i = 0; i < n; i++)
         {
           endo.clear();
           equations[i]->collectEndogenous(endo);
           for (set<pair<int, int> >::const_iterator it = endo.begin(); it != endo.end(); it++)
-            endo_eqs_incidence.insert(make_pair(i, it->first));
+            tmp_normalized_contemporaneous_jacobian[make_pair(i, it->first)] = 1;
         }
-
-      try
+      check = computeNormalization(tmp_normalized_contemporaneous_jacobian, true);
+      if (check)
         {
-          computeNormalization(endo_eqs_incidence);
-        }
-      catch (NormalizationException &e)
-        {
-          cerr << "ERROR: Could not normalize the model even with zero cutoff. Variable "
-               << symbol_table.getName(e.symb_id)
-               << " is not in the maximum cardinality matching." << endl;
-          exit(EXIT_FAILURE);
+          // Update the jacobian matrix
+          for (jacob_map::const_iterator it = tmp_normalized_contemporaneous_jacobian.begin(); it != tmp_normalized_contemporaneous_jacobian.end(); it++)
+            {
+              if (static_jacobian.find(make_pair(it->first.first, it->first.second)) == static_jacobian.end())
+                static_jacobian[make_pair(it->first.first, it->first.second)] = 0;
+              if (dynamic_jacobian.find(make_pair(0, make_pair(it->first.first, it->first.second))) == dynamic_jacobian.end())
+                dynamic_jacobian[make_pair(0, make_pair(it->first.first, it->first.second))] = 0;
+              if (contemporaneous_jacobian.find(make_pair(it->first.first, it->first.second)) == contemporaneous_jacobian.end())
+                contemporaneous_jacobian[make_pair(it->first.first, it->first.second)] = 0;
+            }
         }
     }
+
+  if (!check)
+    {
+      cerr << "No normalization could be computed. Aborting." << endl;
+      exit(EXIT_FAILURE);
+    }
 }
 
 void
diff --git a/preprocessor/ModelTree.hh b/preprocessor/ModelTree.hh
index 7b1d7036a399d5729d7fc147e1e09dd0f5c2c229..4614a12c623c995af13252decdace3e03a4c57ff 100644
--- a/preprocessor/ModelTree.hh
+++ b/preprocessor/ModelTree.hh
@@ -136,24 +136,21 @@ protected:
   //! for each block contains pair< max_lag, max_lead>
   t_lag_lead_vector block_lag_lead;
 
-  //! Exception thrown when normalization fails
-  class NormalizationException
-  {
-  public:
-    //! A variable missing from the maximum cardinal matching
-    int symb_id;
-    NormalizationException(int symb_id_arg) : symb_id(symb_id_arg)
-    {
-    }
-  };
-
   //! Compute the matching between endogenous and variable using the jacobian contemporaneous_jacobian
-  /*! \param endo_eqs_incidence A set indicating which endogenous appear in which equation. First element of pairs is equation number, second is type specific endo ID */
-  void computeNormalization(const set<pair<int, int> > &endo_eqs_incidence) throw (NormalizationException);
+  /*!
+    \param contemporaneous_jacobian Jacobian used as an incidence matrix: all elements declared in the map (even if they are zero), are used as vertices of the incidence matrix
+    \return True if a complete normalization has been achieved
+  */
+  bool computeNormalization(const jacob_map &contemporaneous_jacobian, bool verbose);
+
   //! Try to compute the matching between endogenous and variable using a decreasing cutoff
-  /*! applied to the jacobian contemporaneous_jacobian and stop when a matching is found.*/
-  /*! if no matching is found with a cutoff close to zero an error message is printout */
-  void computePossiblySingularNormalization(const jacob_map &contemporaneous_jacobian, bool try_symbolic);
+  /*!
+    Applied to the jacobian contemporaneous_jacobian and stop when a matching is found.
+    If no matching is found using a strictly positive cutoff, then a zero cutoff is applied (i.e. use a symbolic normalization); in that case, the method adds zeros in the jacobian matrices to reflect all the edges in the symbolic incidence matrix.
+    If no matching is found with a zero cutoff close to zero an error message is printout.
+  */
+  void computeNonSingularNormalization(jacob_map &contemporaneous_jacobian, double cutoff, jacob_map &static_jacobian, dynamic_jacob_map &dynamic_jacobian);
+
   //! Try to normalized each unnormalized equation (matched endogenous variable only on the LHS)
   void computeNormalizedEquations(multimap<int, int> &endo2eqs) const;
   //! Evaluate the jacobian and suppress all the elements below the cutoff
diff --git a/preprocessor/StaticModel.cc b/preprocessor/StaticModel.cc
index 784403e72fdbe767d50d7e67f6e9fd2b14e15b1a..6bce4e2b1a50a7ff8b03515aeb7c728a42957139 100644
--- a/preprocessor/StaticModel.cc
+++ b/preprocessor/StaticModel.cc
@@ -45,7 +45,7 @@ StaticModel::StaticModel(SymbolTable &symbol_table_arg,
 }
 
 void
-StaticModel::compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const
+StaticModel::compileDerivative(ofstream &code_file, int eq, int symb_id, map_idx_type &map_idx) const
 {
   first_derivatives_type::const_iterator it = first_derivatives.find(make_pair(eq, symbol_table.getID(eEndogenous, symb_id)));
   if (it != first_derivatives.end())
@@ -539,7 +539,7 @@ StaticModel::writeModelEquationsCodeOrdered(const string file_name, const string
             {
             case SOLVE_BACKWARD_SIMPLE:
             case SOLVE_FORWARD_SIMPLE:
-              compileDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), 0, map_idx);
+              compileDerivative(code_file, getBlockEquationID(block, 0), getBlockVariableID(block, 0), map_idx);
               {
                 FSTPG_ fstpg(0);
                 fstpg.write(code_file);
@@ -727,7 +727,7 @@ StaticModel::computingPass(const eval_context_type &eval_context, bool no_tmp_te
 
       evaluateAndReduceJacobian(eval_context, contemporaneous_jacobian, static_jacobian, dynamic_jacobian, cutoff, false);
 
-      computePossiblySingularNormalization(contemporaneous_jacobian, cutoff == 0);
+      computeNonSingularNormalization(contemporaneous_jacobian, cutoff, static_jacobian, dynamic_jacobian);
 
       computePrologueAndEpilogue(static_jacobian, equation_reordered, variable_reordered, prologue, epilogue);
 
diff --git a/preprocessor/StaticModel.hh b/preprocessor/StaticModel.hh
index 5b82ab54107d8c02ad812c62677be33ff44c22ca..1dbfb07d133455022a87dfa3f1341b0fcc967c93 100644
--- a/preprocessor/StaticModel.hh
+++ b/preprocessor/StaticModel.hh
@@ -74,7 +74,7 @@ private:
 
   void computeTemporaryTermsOrdered();
   //! Write derivative code of an equation w.r. to a variable
-  void compileDerivative(ofstream &code_file, int eq, int symb_id, int lag, map_idx_type &map_idx) const;
+  void compileDerivative(ofstream &code_file, int eq, int symb_id, map_idx_type &map_idx) const;
   //! Write chain rule derivative code of an equation w.r. to a variable
   void compileChainRuleDerivative(ofstream &code_file, int eq, int var, int lag, map_idx_type &map_idx) const;
 
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 3ead5195123b6f9ba29d46bf5e484c89efee358c..41322bb311deaa0912b29254588d590c62c41472 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -13,6 +13,8 @@ OCTAVE_MODS = \
 	predetermined_variables.mod \
 	block_bytecode/fs2000_simk.mod \
 	block_bytecode/fs2000_lu.mod \
+	block_bytecode/fs2000_bicgstab.mod \
+	block_bytecode/fs2000_optpath.mod \
 	block_bytecode/fs2000_bytecode.mod \
 	block_bytecode/ramst.mod \
 	block_bytecode/ireland.mod \
@@ -48,7 +50,9 @@ MODS = $(OCTAVE_MODS) \
 	AIM/ls2003_2L0L.mod \
 	AIM/ls2003_2L0L_AIM.mod \
 	AIM/ls2003_2L2L.mod \
-	AIM/ls2003_2L2L_AIM.mod
+	AIM/ls2003_2L2L_AIM.mod \
+	block_bytecode/fs2000_gmres.mod \
+	block_bytecode/ramst_a.mod
 
 EXTRA_DIST = $(MODS) \
 	run_test_octave.m \