diff --git a/matlab/dynamic_backward_model_for_simulation.m b/matlab/dynamic_backward_model_for_simulation.m
new file mode 100644
index 000000000..4947aa669
--- /dev/null
+++ b/matlab/dynamic_backward_model_for_simulation.m
@@ -0,0 +1,40 @@
+function [r, J] = dynamic_backward_model_for_simulation(z, dynamicmodel, ylag, x, params, steady_state, it_)
+
+% Copyright (C) 2017 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/>.
+
+% Get indices of the variables appearing at time t.
+% NOTE: It is assumed that all variables appear at time t in the model.
+idy = length(ylag)+(1:length(z));
+
+% Build y vector to be passed to the dynamic model.
+y = zeros(length(ylag)+length(z), 1);
+y(1:length(ylag)) = ylag;
+y(idy) = z;
+
+if nargout>1
+    % Compute residuals and jacobian of the full dynamic model.
+    [r, Jacobian] = feval(dynamicmodel, y, x, params, steady_state, it_);
+else
+    % Compute residuals and return.
+    r = feval(dynamicmodel, y, x, params, steady_state, it_);
+    return
+end
+
+% If the jacobian is computed, remove the columns related to the innovations
+% and the variables appearing at time t-1.
+J = Jacobian(:,idy);
\ No newline at end of file
diff --git a/matlab/simul_backward_nonlinear_model.m b/matlab/simul_backward_nonlinear_model.m
index 1de354e7d..fbd6241d1 100644
--- a/matlab/simul_backward_nonlinear_model.m
+++ b/matlab/simul_backward_nonlinear_model.m
@@ -70,7 +70,6 @@ else
 end
 
 % Get usefull vector of indices.
-ny0 = nnz(DynareModel.lead_lag_incidence(2,:));
 ny1 = nnz(DynareModel.lead_lag_incidence(1,:));
 iy1 = find(DynareModel.lead_lag_incidence(1,:)>0);
 idx = 1:DynareModel.endo_nbr;
@@ -79,6 +78,7 @@ hdx = 1:ny1;
 
 % Get the name of the dynamic model routine.
 model_dynamic = str2func([DynareModel.fname,'_dynamic']);
+model_dynamic_s = str2func('dynamic_backward_model_for_simulation');
 
 % initialization of vector y.
 y = NaN(length(idx)+ny1,1);
@@ -99,14 +99,9 @@ Y = DynareOutput.endo_simul;
 
 % Simulations (call a Newton-like algorithm for each period).
 for it = 2:sample_size+1
-    y(jdx) = Y(:,it-1);                       % A good guess for the initial conditions is the previous values for the endogenous variables.
-    y(hdx) = y(jdx(iy1));                     % Set lagged variables.
-    z = trust_region(model_dynamic, y, idx, jdx, 1, DynareOptions.gstep, ...
-                    DynareOptions.solve_tolf,DynareOptions.solve_tolx, ...
-                    DynareOptions.simul.maxit,DynareOptions.debug, ...
-                    DynareOutput.exo_simul, DynareModel.params, ...
-                    DynareOutput.steady_state, it);
-    Y(:,it) = z(jdx);
+    ylag = Y(iy1,it-1);                   % Set lagged variables.
+    y = Y(:,it-1);                        % A good guess for the initial conditions is the previous values for the endogenous variables.
+    Y(:,it) = dynare_solve(model_dynamic_s, y, DynareOptions, model_dynamic, ylag, DynareOutput.exo_simul, DynareModel.params, DynareOutput.steady_state, it);
 end
 
 DynareOutput.endo_simul = Y;
\ No newline at end of file