diff --git a/dynare++/dynare_simul/dynare_simul.m b/dynare++/dynare_simul/dynare_simul.m
new file mode 100644
index 0000000000000000000000000000000000000000..095d417b7a3489fa2c4e5390209c4d00a9df889c
--- /dev/null
+++ b/dynare++/dynare_simul/dynare_simul.m
@@ -0,0 +1,176 @@
+%
+% SYNOPSIS
+%
+% r = dynare_simul(name, shocks)
+% r = dynare_simul(name, prefix, shocks)
+% r = dynare_simul(name, shocks, start)
+% r = dynare_simul(name, prefix, shocks, start)
+%
+% name name of MAT-file produced by dynare++
+% prefix prefix of variables in the MAT-file
+% shocks matrix of shocks
+% start zero period value
+%
+% Note that this file requires the dynare_simul_ DLL to be in the path.
+% This DLL is distributed with Dynare, under the mex/matlab or mex/octave
+% subdirectory.
+%
+% SEMANTICS
+%
+% The command reads a decision rule from the MAT-file having the given
+% prefix. Then it starts simulating the decision rule with zero time value
+% equal to the given start. It uses the given shocks for the simulation. If
+% the start is not given, the state about which the decision rule is
+% centralized is taken (called fix point, or stochastic steady state, take
+% your pick).
+%
+% prefix Use the prefix with which you called dynare++, the default
+% prefix in dynare++ is 'dyn'.
+% shocks Number of rows must be a number of exogenous shocks,
+% number of columns gives the number of simulated
+% periods. NaNs and Infs in the matrix are substitued by
+% draws from the normal distribution using the covariance
+% matrix given in the model file.
+% start Vector of endogenous variables in the ordering given by
+% <prefix>_vars.
+%
+% Seed for random generator is derived from calling rand(1,1). Therefore,
+% seeding can be controlled with rand('state') and rand('state',some_seed).
+%
+% EXAMPLES
+%
+% All examples suppose that the prefix is 'dyn' and that your_model.mat
+% has been loaded into Matlab.
+%
+% 1. response to permanent negative shock to the third exo var EPS3 for
+% 100 periods
+%
+% shocks = zeros(4,100); % 4 exogenous variables in the model
+% shocks(dyn_i_EPS3,:) = -0.1; % the permanent shock to EPS3
+% r = dynare_simul('your_model.mat',shocks);
+%
+% 2. one stochastic simulation for 100 periods
+%
+% shocks = zeros(4,100)./0; % put NaNs everywhere
+% r = dynare_simul('your_model.mat',shocks);
+%
+% 3. one stochastic simulation starting at 75% undercapitalized economy
+%
+% shocks = zeros(4,100)./0; % put NaNs everywhere
+% ystart = dyn_ss; % get copy of DR fix point
+% ystart(dyn_i_K) = 0.75*dyn_ss(dyn_i_K); % scale down the capital
+% r = dynare_simul('your_model.mat',shocks,ystart);
+%
+%
+% SEE ALSO
+%
+% "DSGE Models with Dynare++. A Tutorial.", Ondra Kamenik, 2005
+
+% Copyright (C) 2005-2011, Ondra Kamenik
+% Copyright (C) 2020, Dynare Team
+
+
+function r = dynare_simul(varargin)
+
+if ~exist('dynare_simul_','file')
+ error('Can''t find dynare_simul_ DLL in the path. The simplest way to add it is to run Dynare once in this session.')
+end
+
+% get the file name and load data
+fname = varargin{1};
+load(fname);
+
+% set prefix, shocks, ystart
+if ischar(varargin{2})
+ prefix = varargin{2};
+ if length(varargin) == 3
+ shocks = varargin{3};
+ ystart = NaN;
+ elseif length(varargin) == 4
+ shocks = varargin{3};
+ ystart = varargin{4};
+ else
+ error('Wrong number of parameters.');
+ end
+else
+ prefix = 'dyn';
+ if length(varargin) == 2
+ shocks = varargin{2};
+ ystart = NaN;
+ elseif length(varargin) == 3
+ shocks = varargin{2};
+ ystart = varargin{3};
+ else
+ error('Wrong number of parameters.');
+ end
+end
+
+% load all needed variables but prefix_g_*
+if exist([prefix '_nstat'],'var')
+ nstat = eval([prefix '_nstat']);
+else
+ error(['Could not find variable ' prefix '_nstat in workspace']);
+end
+if exist([prefix '_npred'],'var')
+ npred = eval([prefix '_npred']);
+else
+ error(['Could not find variable ' prefix '_npred in workspace']);
+end
+if exist([prefix '_nboth'],'var')
+ nboth = eval([prefix '_nboth']);
+else
+ error(['Could not find variable ' prefix '_nboth in workspace']);
+end
+if exist([prefix '_nforw'],'var')
+ nforw = eval([prefix '_nforw']);
+else
+ error(['Could not find variable ' prefix '_nforw in workspace']);
+end
+if exist([prefix '_ss'],'var')
+ ss = eval([prefix '_ss']);
+else
+ error(['Could not find variable ' prefix '_ss in workspace']);
+end
+if exist([prefix '_vcov_exo'],'var')
+ vcov_exo = eval([prefix '_vcov_exo']);
+else
+ error(['Could not find variable ' prefix '_vcov_exo in workspace']);
+end
+nexog = size(vcov_exo,1);
+
+if isnan(ystart)
+ ystart = ss;
+end
+
+% newer version of dynare++ doesn't return prefix_g_0, we make it here if
+% it does not exist in workspace
+g_zero = [prefix '_g_0'];
+if ~exist(g_zero,'var')
+ dr.g_0=zeros(nstat+npred+nboth+nforw,1);
+else
+ dr.g_0=eval(g_zero);
+end
+
+% make derstr a string of comma seperated existing prefix_g_*
+order = 1;
+cont = 1;
+while cont == 1
+ g_ord = [prefix '_g_' num2str(order)];
+ if exist(g_ord,'var')
+ dr.(['g_' num2str(order)])=eval(g_ord);
+ order = order + 1;
+ else
+ cont = 0;
+ end
+end
+
+% set seed
+seed = ceil(10000*rand(1,1));
+
+% call dynare_simul_
+[err,r]=dynare_simul_(order-1,nstat,npred,nboth,nforw,...
+ nexog,ystart,shocks,vcov_exo,seed,ss,dr);
+
+if err
+ error('Simulation failed')
+end
\ No newline at end of file