diff --git a/matlab/dr1.m b/matlab/dr1.m
index f872a6423709bbebddaaf3ce5d461f92c509df94..a8d0c4c1d2c7dc24793a334b513665959d428cbd 100644
--- a/matlab/dr1.m
+++ b/matlab/dr1.m
@@ -79,7 +79,7 @@ if options_.ramsey_policy
     
     k1 = find(M_.orig_model.lead_lag_incidence');
     y = repmat(oo_.dr.ys(1:M_.orig_model.endo_nbr),1,M_.orig_model.maximum_lag+M_.orig_model.maximum_lead+1);
-    [f,fJ,fh] = feval([M_.fname '_dynamic'],y(k1),zeros(1,M_.exo_nbr));
+    [f,fJ,fh] = feval([M_.fname '_dynamic'],y(k1),zeros(1,M_.exo_nbr), M_.params, it_);
 
     % looking for dynamic variables that are both in the original model
     % and in the optimal policy model
@@ -94,7 +94,7 @@ if options_.ramsey_policy
     lead_lag_incidence1 = M_.orig_model.lead_lag_incidence;
     y = repmat(oo_.dr.ys(1:M_.orig_model.endo_nbr),1,M_.orig_model.maximum_lag+M_.orig_model.maximum_lead+1);
     k1 = find(M_.orig_model.lead_lag_incidence');
-    [f,fj,fh] = feval([M_.fname '_dynamic'],y(k1),zeros(1,M_.exo_nbr));
+    [f,fj,fh] = feval([M_.fname '_dynamic'],y(k1),zeros(1,M_.exo_nbr), M_.params, it_);
     nrj = size(fj,1); 
     
     iy = M_.lead_lag_incidence;
@@ -174,11 +174,11 @@ else
   z = z(iyr0) ;
   if options_.order == 1
     [junk,jacobia_] = feval([M_.fname '_dynamic'],z,[oo_.exo_simul ...
-		    oo_.exo_det_simul]);
+		    oo_.exo_det_simul], M_.params, it_);
   elseif options_.order == 2
     [junk,jacobia_,hessian] = feval([M_.fname '_dynamic'],z,...
 				    [oo_.exo_simul ...
-		    oo_.exo_det_simul]);
+		    oo_.exo_det_simul], M_.params, it_);
   end
 end
 
diff --git a/matlab/ff1_.m b/matlab/ff1_.m
index 36dd171bb0dc259ed2585ed1c31a7822b3bd1789..2a645ccf08624a91e7958ef9ade5da4575ebbbc1 100644
--- a/matlab/ff1_.m
+++ b/matlab/ff1_.m
@@ -20,7 +20,7 @@ global it_ M_ oo_
 n1 = size(x,1) - M_.exo_nbr;
 oo_.exo_simul(it_+M_.maximum_lag-M_.maximum_lag,:) = x(n1+1:end)';
 fh = str2func([M_.fname '_static']);
-y=feval(fh,x(1:n1),oo_.exo_simul);
+y=feval(fh,x(1:n1),oo_.exo_simul, M_.params);
 
 
 
diff --git a/matlab/homotopy1.m b/matlab/homotopy1.m
index fb38691b6b63165fe1c8f9fc191790da2b51cb48..fe26bb9b15e1af78e90c1390cb4b7ec87b855caa 100644
--- a/matlab/homotopy1.m
+++ b/matlab/homotopy1.m
@@ -66,7 +66,7 @@ function homotopy1(values, step_nbr)
                                             oo_.steady_state,...
                                             options_.jacobian_flag, ...	    
                                             [oo_.exo_steady_state; ...
-                        oo_.exo_det_steady_state]);
+                        oo_.exo_det_steady_state], M_.params);
   
     if check
       error('HOMOTOPY didn''t succeed')
diff --git a/matlab/homotopy2.m b/matlab/homotopy2.m
index 813d8c1465370f45b082604dc5b21d0475772672..106a62256bc35fea74a787ea24fae26bea2afaef 100755
--- a/matlab/homotopy2.m
+++ b/matlab/homotopy2.m
@@ -80,7 +80,7 @@ function homotopy2(values, step_nbr)
                                               oo_.steady_state,...
                                               options_.jacobian_flag, ...	    
                                               [oo_.exo_steady_state; ...
-                          oo_.exo_det_steady_state]);
+                          oo_.exo_det_steady_state], M_.params);
       
       if check
         error('HOMOTOPY didn''t succeed')
diff --git a/matlab/homotopy3.m b/matlab/homotopy3.m
index 400ab3ec0ab19011f26dddc0c4a6a64d579a0818..d94b55c9c1faa722ec2dfecd14d97177bf4f7c4f 100644
--- a/matlab/homotopy3.m
+++ b/matlab/homotopy3.m
@@ -75,7 +75,7 @@ function homotopy3(values, step_nbr)
                                             oo_.steady_state,...
                                             options_.jacobian_flag, ...	    
                                             [oo_.exo_steady_state; ...
-                        oo_.exo_det_steady_state]);
+                        oo_.exo_det_steady_state], M_.params);
   
     if check
       inc = inc/2;
diff --git a/matlab/initial_estimation_checks.m b/matlab/initial_estimation_checks.m
index 779db48e4eeaba091e88c6f0c3d8ded46229d97b..b29b364328590a49085bb7b85fa4464a7e522259 100644
--- a/matlab/initial_estimation_checks.m
+++ b/matlab/initial_estimation_checks.m
@@ -60,7 +60,7 @@ function initial_estimation_checks(xparam1,gend,data)
 	check1 = max(abs(feval([M_.fname '_static'],...
 			       oo_.steady_state,...
 			       [oo_.exo_steady_state; ...
-		    oo_.exo_det_steady_state]))) > options_.dynatol ;
+              oo_.exo_det_steady_state], M_.params))) > options_.dynatol ;
 	if check1
 	  error(['The seadystate values returned by ' M_.fname ...
 		 '_steadystate.m don''t solve the static model!' ])
diff --git a/matlab/osr1.m b/matlab/osr1.m
index e6230af90ca1327e9255292828b5ddff25a6fc40..58c7d291f207fdf75352e5dbd0557449ab1a2236 100644
--- a/matlab/osr1.m
+++ b/matlab/osr1.m
@@ -26,8 +26,8 @@ function osr1(i_params,i_var,weights)
   
   % check if ys is steady state
   fh = str2func([M_.fname '_static']);
-  if max(abs(feval(fh,oo_.steady_state,exe))) > options_.dynatol
-    [oo_.dr.ys, check] = dynare_solve([M_.fname '_static'],oo_.steady_state,1,exe);
+  if max(abs(feval(fh, oo_.steady_state, exe, M_.params))) > options_.dynatol
+    [oo_.dr.ys, check] = dynare_solve([M_.fname '_static'], oo_.steady_state, 1, exe, M_.params);
     if check
       error('OLR: convergence problem in DYNARE_SOLVE')
     end
diff --git a/matlab/ramsey_dynamic.m b/matlab/ramsey_dynamic.m
index 1635a680d9216f2d6b8d445c96ec1f117f44edef..46029b5cb980cc5d59bf2652b8afdef5243a42b3 100644
--- a/matlab/ramsey_dynamic.m
+++ b/matlab/ramsey_dynamic.m
@@ -54,12 +54,12 @@ function J = ramsey_dynamic(ys,lbar)
   k = find(i_leadlag');
   it_ = 1;
   % retrieving derivatives of the objective function
-  [U,Uy,Uyy] = feval([fname '_objective_static'],ys,zeros(1,exo_nbr));
+  [U,Uy,Uyy] = feval([fname '_objective_static'],ys,zeros(1,exo_nbr), M_.params);
   Uy = Uy';
   Uyy = reshape(Uyy,endo_nbr,endo_nbr);
   
   % retrieving derivatives of original model
-  [f,fJ,fH] = feval([fname '_dynamic'],y(k),zeros(1,exo_nbr));
+  [f,fJ,fH] = feval([fname '_dynamic'],y(k),zeros(1,exo_nbr), M_.params, it_);
   instr_nbr = endo_nbr - size(f,1);
   mult_nbr = endo_nbr-instr_nbr;
 
diff --git a/matlab/ramsey_policy.m b/matlab/ramsey_policy.m
index ec3a4f73796c40f526ce013651b9e78c8af0635f..5eaf17c950eeb2899016257112eeab51863d00d9 100644
--- a/matlab/ramsey_policy.m
+++ b/matlab/ramsey_policy.m
@@ -43,7 +43,7 @@ function info = ramsey_policy(var_list)
   ghs22 = [ghs2(nstatic+(1:npred),:); zeros(nspred-npred,size(ghs2,2))];
   
   fname = [M_.fname '_objective_static'];
-  [ubar,uj,uh] = feval(fname,dr.ys(1:endo_nbr1),zeros(exo_nbr1,1));
+  [ubar,uj,uh] = feval(fname,dr.ys(1:endo_nbr1),zeros(exo_nbr1,1), M_.params);
   
   bet = options_.planner_discount;
   wbar = ubar/(1-bet);
diff --git a/matlab/ramsey_static.m b/matlab/ramsey_static.m
index cd3bb13089cf63ae173955ba6f6c44ae13d87490..8ada7db951788ad841670b84a4ff4d59a391c5a2 100644
--- a/matlab/ramsey_static.m
+++ b/matlab/ramsey_static.m
@@ -38,7 +38,7 @@ function [resids, rJ,mult] = ramsey_static(x)
   
   % value and Jacobian of objective function
   ex = zeros(1,M_.exo_nbr);
-  [U,Uy,Uyy] = feval([fname '_objective_static'],x(i_endo),ex);
+  [U,Uy,Uyy] = feval([fname '_objective_static'],x(i_endo),ex, M_.params);
   Uy = Uy';
   Uyy = reshape(Uyy,endo_nbr,endo_nbr);
   
@@ -47,7 +47,7 @@ function [resids, rJ,mult] = ramsey_static(x)
   k = find(i_lag');
   it_ = 1;
 %  [f,fJ,fH] = feval([fname '_dynamic'],y(k),ex);
-  [f,fJ] = feval([fname '_dynamic'],y(k),ex);
+  [f,fJ] = feval([fname '_dynamic'],y(k),ex, M_.params, it_);
   % indices of Lagrange multipliers
   inst_nbr = endo_nbr - size(f,1);
   i_mult = [endo_nbr+1:2*endo_nbr-inst_nbr]';
diff --git a/matlab/resid.m b/matlab/resid.m
index d3a45d666be4140770f1918308810ba13c2adf19..94e51960e533916b3a11516ea557a1e681ff2e4a 100644
--- a/matlab/resid.m
+++ b/matlab/resid.m
@@ -39,7 +39,7 @@ function resid(period)
   z = zeros(n,period);
   fh = str2func([M_.fname '_dynamic']);
   for it_=M_.maximum_lag+1:period+M_.maximum_lag
-    z(:,it_-M_.maximum_lag) = feval(fh,y(iyr0),oo_.exo_simul);
+    z(:,it_-M_.maximum_lag) = feval(fh,y(iyr0),oo_.exo_simul, M_.params, it_);
     iyr0 = iyr0 + n;
   end
 
diff --git a/matlab/resol.m b/matlab/resol.m
index ee4829a88690533610c066baa9e66bcd77edd103..aaca8f2cec8028732337282873a17e376ea6e078 100644
--- a/matlab/resol.m
+++ b/matlab/resol.m
@@ -60,17 +60,17 @@ if options_.steadystate_flag
 			 [oo_.exo_steady_state; oo_.exo_det_steady_state]);
 else
   % testing if ys isn't a steady state or if we aren't computing Ramsey policy
-  if max(abs(feval(fh,dr.ys,[oo_.exo_steady_state; oo_.exo_det_steady_state]))) ...
+  if max(abs(feval(fh,dr.ys,[oo_.exo_steady_state; oo_.exo_det_steady_state], M_.params))) ...
 	> options_.dynatol & options_.ramsey_policy == 0
     if options_.linear == 0
       % nonlinear models
       [dr.ys,check1] = dynare_solve(fh,dr.ys,options_.jacobian_flag,...
 				    [oo_.exo_steady_state; ...
-		    oo_.exo_det_steady_state]);
+		    oo_.exo_det_steady_state], M_.params);
     else
       % linear models
       [fvec,jacob] = feval(fh,dr.ys,[oo_.exo_steady_state;...
-		    oo_.exo_det_steady_state]);
+		    oo_.exo_det_steady_state], M_.params);
       dr.ys = dr.ys-jacob\fvec;
     end
   end
@@ -78,7 +78,7 @@ end
 % testing for problem
 if check1
   info(1) = 20;
-  resid = feval(fh,ys,oo_.exo_steady_state);
+  resid = feval(fh,ys,oo_.exo_steady_state, M_.params);
   info(2) = resid'*resid; % penalty...
   return
 end
@@ -92,7 +92,7 @@ end
 
 if options_.dr_algo == 1 & options_.order > 1
   dr.ys = dynare_solve('dr2',ys,0,dr);
-  dr.fbias = 2*feval([M_.fname '_static'],dr.ys,oo_.exo_steady_state);
+  dr.fbias = 2*feval([M_.fname '_static'],dr.ys,oo_.exo_steady_state, M_.params);
   [dr, info1] = dr1(dr,check_flag);
   if info1(1)
     info(1) = info(1)+10;
diff --git a/matlab/restricted_steadystate.m b/matlab/restricted_steadystate.m
index e5972e44657968edf6cd6daa264741d651f1d759..7d29374fe433176f5b530785f95f5941a587de7f 100644
--- a/matlab/restricted_steadystate.m
+++ b/matlab/restricted_steadystate.m
@@ -8,7 +8,7 @@ function [sR,sG] = restricted_steadystate(y,x,indx)
   
   ss(indx) = y;
  
-  eval(['[R,G] = ' M_.fname '_static(ss,x);']);
+  eval(['[R,G] = ' M_.fname '_static(ss, x, M_.params);']);
 
   sR = R(inde);
   sG = G(inde,indx);
\ No newline at end of file
diff --git a/matlab/sim1.m b/matlab/sim1.m
index 55d6fac374096293574faf5a3fe557480bedfa74..5e984717d9a7563029f0c2d8d7e6d3ed7720b95d 100644
--- a/matlab/sim1.m
+++ b/matlab/sim1.m
@@ -57,14 +57,14 @@ for iter = 1:options_.maxit
 
 	it_ = it_init ;
 	z = [oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1)] ;
-	[d1,M_.jacobia] = feval([M_.fname '_dynamic'],z,oo_.exo_simul);
+	[d1,M_.jacobia] = feval([M_.fname '_dynamic'],z,oo_.exo_simul, M_.params, it_);
 	M_.jacobia = [M_.jacobia(:,iz) -d1] ;
 	ic = [1:ny] ;
 	icp = iyp ;
 	c (ic,:) = M_.jacobia(:,is)\M_.jacobia(:,isf1) ;
 	for it_ = it_init+(1:options_.periods-1)
 		z = [oo_.endo_simul(iyp,it_-1) ; oo_.endo_simul(:,it_) ; oo_.endo_simul(iyf,it_+1)] ;
-		[d1,M_.jacobia] = feval([M_.fname '_dynamic'],z,oo_.exo_simul);
+		[d1,M_.jacobia] = feval([M_.fname '_dynamic'],z,oo_.exo_simul, M_.params, it_);
 		M_.jacobia = [M_.jacobia(:,iz) -d1] ;
 		M_.jacobia(:,[isf nrs]) = M_.jacobia(:,[isf nrs])-M_.jacobia(:,isp)*c(icp,:) ;
 		ic = ic + ny ;
diff --git a/matlab/simk.m b/matlab/simk.m
index d7200114b832220a045d8bfe6c2d96824abcc486..3631b48f2a82daf2fe0d608549ad650d651fda22 100644
--- a/matlab/simk.m
+++ b/matlab/simk.m
@@ -22,7 +22,7 @@ function simk
 global M_ options_ oo_
 global it_ iyr0 ct_ broyden_
 
-func_name = [M_.fname '_static'];
+%func_name = [M_.fname '_static'];
 nk = M_.maximum_endo_lag + M_.maximum_endo_lead + 1 ;
 ny = size(M_.lead_lag_incidence,2) ;
 icc1 = M_.lead_lag_incidence(nk,:) > 0;
@@ -126,10 +126,10 @@ for iter = 1:options_.maxit
     h3 = clock ;
     if broyden_ & iter > 1
       %d1_ = -feval(fh,oo_.endo_simul(iyr));
-      d1 = -feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),z,oo_.exo_simul);
+      d1 = -feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),z,oo_.exo_simul, M_.params, it_);
     else
       %jacob(func_name,oo_.endo_simul(iyr)) ;
-      [d1,M_.jacobia] = feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),oo_.exo_simul);
+      [d1,M_.jacobia] = feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),oo_.exo_simul, M_.params, it_);
       d1 = -d1 ;
     end
     err_f = max(err_f,max(abs(d1)));
@@ -219,10 +219,10 @@ for iter = 1:options_.maxit
   while it_ <= options_.periods+M_.maximum_lag
     if broyden_
       %d1_ = -feval(fh,oo_.endo_simul(iyr));
-       d1 = -feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),z,oo_.exo_simul);
+       d1 = -feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),z,oo_.exo_simul, M_.params, it_);
     else
       %jacob(func_name,oo_.endo_simul(iyr)) ;
-      [d1,M_.jacobia] = feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),oo_.exo_simul);
+      [d1,M_.jacobia] = feval([M_.fname '_dynamic'],oo_.endo_simul(iyr),oo_.exo_simul, M_.params, it_);
       d1 = -d1 ;
     end
     err_f = max(err_f,max(abs(d1)));
diff --git a/matlab/steady_.m b/matlab/steady_.m
index 08cba218ab9e4905e75f73ceaf99bae214978635..34499332b59b317013f8a59111a8490a5b052510 100644
--- a/matlab/steady_.m
+++ b/matlab/steady_.m
@@ -32,7 +32,7 @@ function steady_()
 	check1 = max(abs(feval([M_.fname '_static'],...
 			       oo_.steady_state,...
 			       [oo_.exo_steady_state; ...
-		    oo_.exo_det_steady_state]))) > options_.dynatol ;
+		    oo_.exo_det_steady_state], M_.params))) > options_.dynatol ;
 	if check1
 	  error(['The seadystate values returned by ' M_.fname ...
 		'_steadystate.m don''t solve the static model!' ])
@@ -56,7 +56,7 @@ function steady_()
 				     oo_.steady_state,...
 				     options_.jacobian_flag, ...	    
 			             [oo_.exo_steady_state; ...
-		                      oo_.exo_det_steady_state]);
+                    oo_.exo_det_steady_state], M_.params);
   end
 
   if check ~= 0
diff --git a/preprocessor/ModelTree.cc b/preprocessor/ModelTree.cc
index 8c76a171be4b11992e910b7d742f851cab4c6e26..8090ef0ed0f7dfa5708264956e78bd784af43a18 100644
--- a/preprocessor/ModelTree.cc
+++ b/preprocessor/ModelTree.cc
@@ -1530,7 +1530,7 @@ ModelTree::writeStaticMFile(const string &static_basename) const
       exit(-1);
     }
   // Writing comments and function definition command
-  mStaticModelFile << "function [residual, g1, g2] = " << static_basename << "( y, x )\n";
+  mStaticModelFile << "function [residual, g1, g2] = " << static_basename << "(y, x, params)" << endl;
   mStaticModelFile << interfaces::comment()+"\n"+interfaces::comment();
   mStaticModelFile << "Status : Computes static model for Dynare\n" << interfaces::comment() << "\n";
   mStaticModelFile << interfaces::comment();
@@ -1557,7 +1557,7 @@ ModelTree::writeDynamicMFile(const string &dynamic_basename) const
       cerr << "Error: Can't open file " << filename << " for writing" << endl;
       exit(-1);
     }
-  mDynamicModelFile << "function [residual, g1, g2, g3] = " << dynamic_basename << "(y, x)\n";
+  mDynamicModelFile << "function [residual, g1, g2, g3] = " << dynamic_basename << "(y, x, params, it_)" << endl;
   mDynamicModelFile << interfaces::comment()+"\n"+interfaces::comment();
   mDynamicModelFile << "Status : Computes dynamic model for Dynare\n" << interfaces::comment() << "\n";
   mDynamicModelFile << interfaces::comment();
@@ -1590,9 +1590,7 @@ ModelTree::writeStaticCFile(const string &static_basename) const
                    << endl
                    << " */" << endl
                    << "#include <math.h>" << endl
-                   << "#include \"mex.h\"" << endl
-    // A global variable for model parameters
-                   << "double *params;" << endl;
+                   << "#include \"mex.h\"" << endl;
 
   // Writing the function Static
   writeStaticModel(mStaticModelFile);
@@ -1601,9 +1599,8 @@ ModelTree::writeStaticCFile(const string &static_basename) const
   mStaticModelFile << "/* The gateway routine */" << endl
                    << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
                    << "{" << endl
-                   << "  double *y, *x;" << endl
+                   << "  double *y, *x, *params;" << endl
                    << "  double *residual, *g1;" << endl
-                   << "  mxArray *M_;" << endl
                    << endl
                    << "  /* Create a pointer to the input matrix y. */" << endl
                    << "  y = mxGetPr(prhs[0]);" << endl
@@ -1611,6 +1608,9 @@ ModelTree::writeStaticCFile(const string &static_basename) const
                    << "  /* Create a pointer to the input matrix x. */" << endl
                    << "  x = mxGetPr(prhs[1]);" << endl
                    << endl
+                   << "  /* Create a pointer to the input matrix params. */" << endl
+                   << "  params = mxGetPr(prhs[2]);" << endl
+                   << endl
                    << "  residual = NULL;" << endl
                    << "  if (nlhs >= 1)" << endl
                    << "  {" << endl
@@ -1629,15 +1629,8 @@ ModelTree::writeStaticCFile(const string &static_basename) const
                    << "      g1 = mxGetPr(plhs[1]);" << endl
                    << "  }" << endl
                    << endl
-                   << "  /* Gets model parameters from global workspace of Matlab */" << endl
-                   << "  M_ = mexGetVariable(\"global\",\"M_\");" << endl
-                   << "  if (M_ == NULL ){" << endl
-                   << "     mexPrintf(\"Global variable not found : \");" << endl
-                   << "     mexErrMsgTxt(\"M_ \\n\");" << endl
-                   << "  }" << endl
-                   << "  params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"params\")));" << endl
                    << "  /* Call the C Static. */" << endl
-                   << "  Static(y, x, residual, g1);" << endl
+                   << "  Static(y, x, params, residual, g1);" << endl
                    << "}" << endl;
 
   mStaticModelFile.close();
@@ -1663,12 +1656,7 @@ ModelTree::writeDynamicCFile(const string &dynamic_basename) const
                     << endl
                     << " */" << endl
                     << "#include <math.h>" << endl
-                    << "#include \"mex.h\"" << endl
-    // A global variable for model parameters
-                    << "double *params;" << endl
-  // A global variable for it_
-                    << "int it_;" << endl
-                    << "int nb_row_x;" << endl;
+                    << "#include \"mex.h\"" << endl;
 
   // Writing the function body
   writeDynamicModel(mDynamicModelFile);
@@ -1677,15 +1665,22 @@ ModelTree::writeDynamicCFile(const string &dynamic_basename) const
   mDynamicModelFile << "/* The gateway routine */" << endl
                     << "void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])" << endl
                     << "{" << endl
-                    << "  double *y, *x;" << endl
+                    << "  double *y, *x, *params;" << endl
                     << "  double *residual, *g1, *g2;" << endl
-                    << "  mxArray *M_;" << endl
+                    << "  int nb_row_x, it_;" << endl
                     << endl
                     << "  /* Create a pointer to the input matrix y. */" << endl
                     << "  y = mxGetPr(prhs[0]);" << endl
                     << endl
                     << "  /* Create a pointer to the input matrix x. */" << endl
                     << "  x = mxGetPr(prhs[1]);" << endl
+                    << endl
+                    << "  /* Create a pointer to the input matrix params. */" << endl
+                    << "  params = mxGetPr(prhs[2]);" << endl
+                    << endl
+                    << "  /* Fetch time index */" << endl
+                    << "  it_ = (int) mxGetScalar(prhs[3]) - 1;" << endl
+                    << endl
                     << "  /* Gets number of rows of matrix x. */" << endl
                     << "  nb_row_x = mxGetM(prhs[1]);" << endl
                     << endl
@@ -1721,18 +1716,8 @@ ModelTree::writeDynamicCFile(const string &dynamic_basename) const
                     << "     g2 = mxGetPr(plhs[2]);" << endl
                     << "  }" << endl
                     << endl
-                    << "  /* Gets model parameters from global workspace of Matlab */" << endl
-                    << "  M_ = mexGetVariable(\"global\",\"M_\");" << endl
-                    << "  if (M_ == NULL)" << endl
-                    << "  {" << endl
-                    << "      mexPrintf(\"Global variable not found : \");" << endl
-                    << "      mexErrMsgTxt(\"M_ \\n\");" << endl
-                    << "  }" << endl
-                    << "  params = mxGetPr(mxGetFieldByNumber(M_, 0, mxGetFieldNumber(M_,\"params\")));" << endl
-                    << "  /* Gets it_ from global workspace of Matlab */" << endl
-                    << "  it_ = (int) mxGetScalar(mexGetVariable(\"global\", \"it_\"))-1;" << endl
                     << "  /* Call the C subroutines. */" << endl
-                    << "  Dynamic(y, x, residual, g1, g2);" << endl
+                    << "  Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2);" << endl
                     << "}" << endl;
   mDynamicModelFile.close();
 }
@@ -1814,9 +1799,6 @@ ModelTree::writeStaticModel(ostream &StaticOutput) const
   // Writing ouputs
   if (mode != eDLLMode)
     {
-      StaticOutput << "global M_ \n";
-      StaticOutput << "if M_.param_nbr > 0\n  params = M_.params;\nend\n";
-
       StaticOutput << "  residual = zeros( " << equations.size() << ", 1);\n";
       StaticOutput << "\n\t"+interfaces::comment()+"\n\t"+interfaces::comment();
       StaticOutput << "Model equations\n\t";
@@ -1854,7 +1836,7 @@ ModelTree::writeStaticModel(ostream &StaticOutput) const
     }
   else
     {
-      StaticOutput << "void Static(double *y, double *x, double *residual, double *g1)" << endl
+      StaticOutput << "void Static(double *y, double *x, double *params, double *residual, double *g1)" << endl
                    << "{" << endl
                    << "  double lhs, rhs;" << endl
         // Writing residual equations
@@ -3419,10 +3401,8 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const
 
   if (mode == eStandardMode)
     {
-      DynamicOutput << "global M_ it_\n";
-      DynamicOutput << "if M_.param_nbr > 0\n  params =  M_.params;\nend\n";
-      DynamicOutput << "\n\t"+interfaces::comment()+"\n\t"+interfaces::comment();
-      DynamicOutput << "Model equations\n\t";
+      DynamicOutput << interfaces::comment() << endl << interfaces::comment();
+      DynamicOutput << "Model equations" << endl;
       DynamicOutput << interfaces::comment() + "\n\n";
       DynamicOutput << "residual = zeros(" << nrows << ", 1);\n";
 
@@ -3466,7 +3446,7 @@ ModelTree::writeDynamicModel(ostream &DynamicOutput) const
     }
   else
     {
-      DynamicOutput << "void Dynamic(double *y, double *x, double *residual, double *g1, double *g2)" << endl
+      DynamicOutput << "void Dynamic(double *y, double *x, int nb_row_x, double *params, int it_, double *residual, double *g1, double *g2)" << endl
                     << "{" << endl
                     << "  double lhs, rhs;" << endl
                     << endl