diff --git a/matlab/dynare_config.m b/matlab/dynare_config.m
index 57f91611342fc0fb5e6fe66024e0f5ea3379fb38..cf671bd6306f1690b69a8bfecaeada75a9200b9e 100644
--- a/matlab/dynare_config.m
+++ b/matlab/dynare_config.m
@@ -86,6 +86,11 @@ if ~exist('OCTAVE_VERSION') && matlab_ver_less_than('7.4')
     addpath([dynareroot '/missing/bsxfun'])
 end
 
+% ilu is missing in old versions of MATLAB and in Octave
+if exist('OCTAVE_VERSION') || matlab_ver_less_than('7.4')
+    addpath([dynareroot '/missing/ilu'])
+end
+
 % nanmean is in Octave Forge Statistics package and in MATLAB Statistics
 % toolbox
 if (exist('OCTAVE_VERSION') && ~user_has_octave_forge_package('statistics')) ...
diff --git a/matlab/solve_one_boundary.m b/matlab/solve_one_boundary.m
index 8ae660889da51e0721be20eb46eb58a0b5611e8f..3e882557571ca4d677a8940c3959a7f335bffd7c 100644
--- a/matlab/solve_one_boundary.m
+++ b/matlab/solve_one_boundary.m
@@ -55,7 +55,7 @@ function [y, info] = solve_one_boundary(fname, y, x, params, steady_state, ...
 %   none.
 %  
 
-% Copyright (C) 1996-2012 Dynare Team
+% Copyright (C) 1996-2013 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -77,7 +77,7 @@ Blck_size=size(y_index_eq,2);
 g2 = [];
 g3 = [];
 correcting_factor=0.01;
-luinc_tol=1e-10;
+ilu_setup.droptol=1e-10;
 max_resa=1e100;
 reduced = 0;
 if(forward_backward)
@@ -316,7 +316,7 @@ for it_=start:incr:finish
                     disp('steady: GMRES ');
                 end
                 while(flag1>0)
-                    [L1, U1]=luinc(g1,luinc_tol);
+                    [L1, U1]=ilu(g1,ilu_setup);
                     [dx,flag1] = gmres(g1,-r,Blck_size,1e-6,Blck_size,L1,U1);
                     if (flag1>0 || reduced)
                         if(flag1==1)
@@ -326,7 +326,7 @@ for it_=start:incr:finish
                         elseif(flag1==3)
                             disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')]);
                         end;
-                        luinc_tol = luinc_tol/10;
+                        ilu_setup.droptol = ilu_setup.droptol/10;
                         reduced = 0;
                     else
                         ya = ya + lambda*dx;
@@ -343,7 +343,7 @@ for it_=start:incr:finish
                     disp('steady: BiCGStab');
                 end
                 while(flag1>0)
-                    [L1, U1]=luinc(g1,luinc_tol);
+                    [L1, U1]=ilu(g1,ilu_setup);
                     phat = ya - U1 \ (L1 \ r);
                     if(is_dynamic)
                         y(it_,y_index_eq) = phat;
@@ -370,7 +370,7 @@ for it_=start:incr:finish
                         elseif(flag1==3)
                             disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Num,'%3d')]);
                         end;
-                        luinc_tol = luinc_tol/10;
+                        ilu_setup.droptol = ilu_setup.droptol/10;
                         reduced = 0;
                     else
                         ya = ya + lambda*dx;
diff --git a/matlab/solve_two_boundaries.m b/matlab/solve_two_boundaries.m
index cde94b9dda7d3a12821514fd5204765432b6ca03..19b1cc392e10ab1c455865c3df31fac38021e3ff 100644
--- a/matlab/solve_two_boundaries.m
+++ b/matlab/solve_two_boundaries.m
@@ -43,7 +43,7 @@ function y = solve_two_boundaries(fname, y, x, params, steady_state, y_index, nz
 %   none.
 %  
 
-% Copyright (C) 1996-2011 Dynare Team
+% Copyright (C) 1996-2013 Dynare Team
 %
 % This file is part of Dynare.
 %
@@ -67,7 +67,7 @@ g2 = [];
 g3 = [];
 Blck_size=size(y_index,2);
 correcting_factor=0.01;
-luinc_tol=1e-10;
+ilu_setup.droptol=1e-10;
 max_resa=1e100;
 Jacobian_Size=Blck_size*(y_kmin+y_kmax_l +periods);
 g1=spalloc( Blck_size*periods, Jacobian_Size, nze*periods);
@@ -255,7 +255,7 @@ while ~(cvg==1 || iter>maxit_),
         elseif(stack_solve_algo==2),
             flag1=1;
             while(flag1>0)
-                [L1, U1]=luinc(g1a,luinc_tol);
+                [L1, U1]=ilu(g1a,ilu_setup);
                 [za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);
                 if (flag1>0 || reduced)
                     if(flag1==1)
@@ -265,7 +265,7 @@ while ~(cvg==1 || iter>maxit_),
                     elseif(flag1==3)
                         disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Size,'%3d')]);
                     end;
-                    luinc_tol = luinc_tol/10;
+                    ilu_setup.droptol = ilu_setup.droptol/10;
                     reduced = 0;
                 else
                     dx = za - ya;
@@ -276,7 +276,7 @@ while ~(cvg==1 || iter>maxit_),
         elseif(stack_solve_algo==3),
             flag1=1;
             while(flag1>0)
-                [L1, U1]=luinc(g1a,luinc_tol);
+                [L1, U1]=ilu(g1a,ilu_setup);
                 [za,flag1] = bicgstab(g1a,b,1e-7,Blck_size*periods,L1,U1);
                 if (flag1>0 || reduced)
                     if(flag1==1)
@@ -286,7 +286,7 @@ while ~(cvg==1 || iter>maxit_),
                     elseif(flag1==3)
                         disp(['Error in simul: GMRES stagnated (Two consecutive iterates were the same.), in block' num2str(Block_Size,'%3d')]);
                     end;
-                    luinc_tol = luinc_tol/10;
+                    ilu_setup.droptol = ilu_setup.droptol/10;
                     reduced = 0;
                 else
                     dx = za - ya;
diff --git a/mex/sources/bytecode/SparseMatrix.cc b/mex/sources/bytecode/SparseMatrix.cc
index c9870603c142563bb6ec5419425d367272126e23..ded9eb64ee2b1db56e279b5bd46a8166f8cdfd36 100644
--- a/mex/sources/bytecode/SparseMatrix.cc
+++ b/mex/sources/bytecode/SparseMatrix.cc
@@ -1,5 +1,5 @@
 /*
- * Copyright (C) 2007-2012 Dynare Team
+ * Copyright (C) 2007-2013 Dynare Team
  *
  * This file is part of Dynare.
  *
@@ -1954,11 +1954,15 @@ SparseMatrix::Solve_Matlab_GMRES(mxArray *A_m, mxArray *b_m, int Size, double sl
   throw FatalExceptionHandling(tmp.str());
 #endif
   int n = mxGetM(A_m);
+  const char *field_names[] = {"droptol"};
+  mwSize dims[1] = { 1 };
+  mxArray *Setup = mxCreateStructArray(1, dims, 1, field_names);
+  mxSetFieldByNumber(Setup, 0, 0, mxCreateDoubleScalar(lu_inc_tol));
   mxArray *lhs0[2];
   mxArray *rhs0[2];
   rhs0[0] = A_m;
-  rhs0[1] = mxCreateDoubleScalar(lu_inc_tol);
-  mexCallMATLAB(2, lhs0, 2, rhs0, "luinc");
+  rhs0[1] = Setup;
+  mexCallMATLAB(2, lhs0, 2, rhs0, "ilu");
   mxArray *L1 = lhs0[0];
   mxArray *U1 = lhs0[1];
   /*[za,flag1] = gmres(g1a,b,Blck_size,1e-6,Blck_size*periods,L1,U1);*/
@@ -2033,11 +2037,15 @@ SparseMatrix::Solve_Matlab_BiCGStab(mxArray *A_m, mxArray *b_m, int Size, double
 {
   unsigned int n = mxGetM(A_m);
   /*[L1, U1]=luinc(g1a,luinc_tol);*/
+  const char *field_names[] = {"droptol"};
+  mwSize dims[1] = { 1 };
+  mxArray *Setup = mxCreateStructArray(1, dims, 1, field_names);
+  mxSetFieldByNumber(Setup, 0, 0, mxCreateDoubleScalar(lu_inc_tol));
   mxArray *lhs0[2];
   mxArray *rhs0[2];
   rhs0[0] = A_m;
-  rhs0[1] = mxCreateDoubleScalar(lu_inc_tol);
-  mexCallMATLAB(2, lhs0, 2, rhs0, "luinc");
+  rhs0[1] = Setup;
+  mexCallMATLAB(2, lhs0, 2, rhs0, "ilu");
   mxArray *L1 = lhs0[0];
   mxArray *U1 = lhs0[1];
   double flags = 2;