diff --git a/mex/sources/k_order_perturbation/tests/first_order.m b/mex/sources/k_order_perturbation/tests/first_order.m
deleted file mode 100644
index c1a2599bd..000000000
--- a/mex/sources/k_order_perturbation/tests/first_order.m
+++ /dev/null
@@ -1,129 +0,0 @@
-function [gy]=first_order(M_, dr, jacobia)
-% Emulation of Dynare++ c++ first_order.cpp for testing pruposes
-
-% Copyright © 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 .
-
-% fd = jacobia_
-% reorder jacobia_
-[fd]=k_reOrderedJacobia(M_, jacobia)
-
-%ypart=dr;
-ypart.ny=M_.endo_nbr;
-ypart.nyss=dr.nfwrd+dr.nboth;
-ypart.nys=dr.npred;
-ypart.npred=dr.npred-dr.nboth;
-ypart.nboth=dr.nboth;
-ypart.nforw=dr.nfwrd;
-ypart.nstat =dr.nstatic
-nu=M_.exo_nbr
-
-off= 1; % = 0 in C
-fyplus = fd(:,off:off+ypart.nyss-1);
-off= off+ypart.nyss;
-fyszero= fd(:,off:off+ypart.nstat-1);
-off= off+ypart.nstat;
-fypzero= fd(:,off:off+ypart.npred-1);
-off= off+ypart.npred;
-fybzero= fd(:,off:off+ypart.nboth-1);
-off= off+ypart.nboth;
-fyfzero= fd(:,off:off+ypart.nforw-1);
-off= off+ypart.nforw;
-fymins= fd(:,off:off+ypart.nys-1);
-off= off+ypart.nys;
-fuzero= fd(:,off:off+nu-1);
-off=off+ nu;
-
-n= ypart.ny+ypart.nboth;
-%TwoDMatrix
-matD=zeros(n,n);
-% matD.place(fypzero,0,0);
-matD(1:n-ypart.nboth,1:ypart.npred)= fypzero;
-% matD.place(fybzero,0,ypart.npred);
-matD(1:n-ypart.nboth,ypart.npred+1:ypart.npred+ypart.nboth)=fybzero;
-% matD.place(fyplus,0,ypart.nys()+ypart.nstat);
-matD(1:n-ypart.nboth,ypart.nys+ypart.nstat+1:ypart.nys+ypart.nstat+ypart.nyss)=fyplus;
-for i=1:ypart.nboth
- matD(ypart.ny()+i,ypart.npred+i)= 1.0;
-end
-
-matE=[fymins, fyszero, zeros(n-ypart.nboth,ypart.nboth), fyfzero; zeros(ypart.nboth,n)];
-% matE.place(fymins;
-% matE.place(fyszero,0,ypart.nys());
-% matE.place(fyfzero,0,ypart.nys()+ypart.nstat+ypart.nboth);
-
-for i= 1:ypart.nboth
- matE(ypart.ny()+i,ypart.nys()+ypart.nstat+i)= -1.0;
-end
-matE=-matE; %matE.mult(-1.0);
-
-% vsl=zeros(n,n);
-% vsr=zeros(n,n);
-% lwork= 100*n+16;
-% work=zeros(1,lwork);
-% bwork=zeros(1,n);
-%int info;
-
-% LAPACK_dgges("N","V","S",order_eigs,&n,matE.getData().base(),&n,
-% matD.getData().base(),&n,&sdim,alphar.base(),alphai.base(),
-% beta.base(),vsl.getData().base(),&n,vsr.getData().base(),&n,
-% work.base(),&lwork,&(bwork[0]),&info);
-
-[matE1,matD1,vsr,sdim,dr.eigval,info] = mjdgges(matE,matD,1);
-
-bk_cond= (sdim==ypart.nys);
-
-% ConstGeneralMatrix z11(vsr,0,0,ypart.nys(),ypart.nys());
-z11=vsr(1:ypart.nys,1:ypart.nys);
-% ConstGeneralMatrix z12(vsr,0,ypart.nys(),ypart.nys(),n-ypart.nys());
-z12=vsr(1:ypart.nys(),ypart.nys+1:end);%, n-ypart.nys);
- % ConstGeneralMatrix z21(vsr,ypart.nys(),0,n-ypart.nys(),ypart.nys());
-z21=vsr(ypart.nys+1:end,1:ypart.nys);
-% ConstGeneralMatrix z22(vsr,ypart.nys(),ypart.nys(),n-ypart.nys(),n-ypart.nys());
-z22=vsr(ypart.nys+1:end,ypart.nys+1:end);
-
-% GeneralMatrix sfder(z12,"transpose");
-sfder=z12';%,"transpose");
- % z22.multInvLeftTrans(sfder);
-sfder=z22'\sfder;
-sfder=-sfder;% .mult(-1);
-
-%s11(matE,0,0,ypart.nys(),ypart.nys());
-s11=matE1(1:ypart.nys,1:ypart.nys);
-% t11=(matD1,0,0,ypart.nys(),ypart.nys());
-t11=matD1(1:ypart.nys,1:ypart.nys);
-dumm=(s11');%,"transpose");
- %z11.multInvLeftTrans(dumm);
-dumm=z11'\dumm;
-preder=(dumm');%,"transpose");
- %t11.multInvLeft(preder);
-preder=t11\preder;
-%preder.multLeft(z11);
-preder= z11*preder;
-
-% gy.place(preder,ypart.nstat,0);
-% gy=(zeros(ypart.nstat,size(preder,2)) ;preder);
-% sder(sfder,0,0,ypart.nstat,ypart.nys());
-sder=sfder(1:ypart.nstat,1:ypart.nys);
-% gy.place(sder,0,0);
-% gy(1:ypart.nstat, 1:ypart.nys)=sder;
-% gy=[sder;preder];
-% fder(sfder,ypart.nstat+ypart.nboth,0,ypart.nforw,ypart.nys());
-fder=sfder(ypart.nstat+ypart.nboth+1:ypart.nstat+ypart.nboth+ypart.nforw,1:ypart.nys);
-% gy.place(fder,ypart.nstat+ypart.nys(),0);
-% gy(ypart.nstat+ypart.nys,1)=fder;
-gy=[sder;preder;fder];
diff --git a/mex/sources/k_order_perturbation/tests/fs2000k.mod b/mex/sources/k_order_perturbation/tests/fs2000k.mod
deleted file mode 100644
index 850b1f938..000000000
--- a/mex/sources/k_order_perturbation/tests/fs2000k.mod
+++ /dev/null
@@ -1,119 +0,0 @@
-// This file replicates the estimation of the CIA model from
-// Frank Schorfheide (2000) "Loss function-based evaluation of DSGE models"
-// Journal of Applied Econometrics, 15, 645-670.
-// the data are the ones provided on Schorfheide's web site with the programs.
-// http://www.econ.upenn.edu/~schorf/programs/dsgesel.ZIP
-// You need to have fsdat.m in the same directory as this file.
-// This file replicates:
-// -the posterior mode as computed by Frank's Gauss programs
-// -the parameter mean posterior estimates reported in the paper
-// -the model probability (harmonic mean) reported in the paper
-// This file was tested with dyn_mat_test_0218.zip
-// the smooth shocks are probably stil buggy
-//
-// The equations are taken from J. Nason and T. Cogley (1994)
-// "Testing the implications of long-run neutrality for monetary business
-// cycle models" Journal of Applied Econometrics, 9, S37-S70.
-// Note that there is an initial minus sign missing in equation (A1), p. S63.
-//
-// Michel Juillard, February 2004
-// Modified for testing k_order_perturbation by GP, Jan-Feb 09
-
-options_.usePartInfo=0;
-options_.use_k_order=0;
-
-//var m m_1 P P_1 c e W R k d n l gy_obs gp_obs Y_obs P_obs y dA P2 c2;
-var m m_1 P P_1 c e W R k d n l gy_obs gp_obs y dA P2 c2;
-varexo e_a e_m;
-
-parameters alp bet gam mst rho psi del;
-
-alp = 0.33;
-bet = 0.99;
-gam = 0.003;
-mst = 1.011;
-rho = 0.7;
-psi = 0.787;
-del = 0.02;
-
-model (use_dll);
-dA = exp(gam+e_a);
-log(m) = (1-rho)*log(mst) + rho*log(m_1(-1))+e_m;
--P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c2(+1)*P2(+1)*m(+1))=0;
-W = l/n;
--(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
-R = P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(-alp)/W;
-1/(c*P)-bet*P*(1-alp)*exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)/(m*l*c(+1)*P(+1)) = 0;
-c+k = exp(-alp*(gam+e_a))*k(-1)^alp*n^(1-alp)+(1-del)*exp(-(gam+e_a))*k(-1);
-P*c = m;
-m-1+d = l;
-e = exp(e_a);
-y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
-gy_obs = dA*y/y(-1);
-gp_obs = (P/P_1(-1))*m_1(-1)/dA;
-//Y_obs/Y_obs(-1) = gy_obs;
-//P_obs/P_obs(-1) = gp_obs;
-P2 = P(+1);
-c2 = c(+1);
-m_1 = m;
-P_1 = P;
-end;
-
-initval;
-m = mst;
-m_1=mst;
-P = 2.25;
-P_1 = 2.25;
-c = 0.45;
-e = 1;
-W = 4;
-R = 1.02;
-k = 6;
-d = 0.85;
-n = 0.19;
-l = 0.86;
-y = 0.6;
-gy_obs = exp(gam);
-gp_obs = exp(-gam);
-dA = exp(gam);
- P2=P;
- c2=c;
-end;
-
-shocks;
-var e_a; stderr 0.014;
-var e_m; stderr 0.005;
-end;
-
-//unit_root_vars P_obs Y_obs;
-
-steady(solve_algo = 2);
-
-check;
-
-estimated_params;
-alp, beta_pdf, 0.356, 0.02;
-bet, beta_pdf, 0.993, 0.002;
-gam, normal_pdf, 0.0085, 0.003;
-mst, normal_pdf, 1.0002, 0.007;
-rho, beta_pdf, 0.129, 0.223;
-psi, beta_pdf, 0.65, 0.05;
-del, beta_pdf, 0.01, 0.005;
-stderr e_a, inv_gamma_pdf, 0.035449, inf;
-stderr e_m, inv_gamma_pdf, 0.008862, inf;
-end;
-
-//varobs P_obs Y_obs;
-varobs gp_obs gy_obs;
-
-steady(solve_algo = 2);
-
-//observation_trends;
-//P_obs (log(mst)-gam);
-//Y_obs (gam);
-//end;
-
-//options_.useAIM = 1;
-estimation(datafile=fsdat,nobs=192,loglinear,mh_replic=2000,
- mode_compute=4,mh_nblocks=2,mh_drop=0.45,mh_jscale=0.65);
-
diff --git a/mex/sources/k_order_perturbation/tests/k_order_test_main.cc b/mex/sources/k_order_perturbation/tests/k_order_test_main.cc
deleted file mode 100644
index 2c8612f79..000000000
--- a/mex/sources/k_order_perturbation/tests/k_order_test_main.cc
+++ /dev/null
@@ -1,328 +0,0 @@
-/*
- * Copyright © 2008-2011 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 .
- */
-
-/*************************************
- * This main() is for testing k_order DLL entry point by linking to
- * the k_ord library statically and passing its hard-coded data:
- * parameters, covar, ysteady and the variable names from fs2000a.mod model
- * The main has been derived from mxFunction used for K-Order DLL
- ***************************************/
-
-// #include "stdafx.h"
-#include "dynamic_dll.h"
-#include "k_ord_dynare.h"
-
-int
-main(int argc, char* argv[])
-{
-
- double qz_criterium = 1 + 1e-6;
- const int check_flag = 0;
- const char* fName = "./fs2000k"; // mxArrayToString(mFname);
- const char* dfExt = ".mexa64"; // Dynamic file extension, e.g.".dll";
-
-#ifdef DEBUG
- mexPrintf("k_order_perturbation: check_flag = %d , fName = %s .\n", check_flag, fName);
-#endif
- int kOrder = 2;
- int npar = 7; //(int)mxGetM(mxFldp);
- double dparams[7] = {0.3300, 0.9900, 0.0030, 1.0110, 0.7000, 0.7870, 0.0200};
- Vector* modParams = new Vector(dparams, npar);
-
-#ifdef DEBUG
- mexPrintf("k_ord_perturbation: nParams=%d .\n", npar);
- for (int i = 0; i < npar; i++)
- {
- mexPrintf("k_ord_perturbation: dParams[%d]= %g.\n", i, dparams + i * (sizeof(double)));
- }
- for (int i = 0; i < npar; i++)
- {
- mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]);
- }
-
-#endif
-
- double d2Dparams[4] = {//(double *) mxGetData(mxFldp);
- 0.1960e-3, 0.0, 0.0, 0.0250e-3};
- npar = 2; //(int)mxGetN(mxFldp);
- TwoDMatrix* vCov = new TwoDMatrix(npar, npar, (d2Dparams));
- double dYSparams[16] = {
- // 27 mxGetData(mxFldp);
- // 1.0110, 2.2582, 5.8012, 0.5808,
- 1.0110, 2.2582, 0.4477, 1.0000, 4.5959, 1.0212, 5.8012, 0.8494, 0.1872,
- 0.8604, 1.0030, 1.0080, 0.5808, 1.0030, 2.2582, 0.4477
- //, 1.0110, 2.2582, 0.4477, 1.0000, 0.1872, 2.2582, 0.4477
- };
- const int nSteady = 16; // 27 //31;//29, 16 (int)mxGetM(mxFldp);
- Vector* ySteady = new Vector(dYSparams, nSteady);
-
- double nnzd[3] = {77, 217, 0};
- const Vector* NNZD = new Vector(nnzd, 3);
-
- // mxFldp = mxGetField(dr, 0,"nstatic" );
- const int nStat = 7; //(int)mxGetScalar(mxFldp);
- // mxFldp = mxGetField(dr, 0,"npred" );
- const int nPred = 2; // 6 - nBoth (int)mxGetScalar(mxFldp);
- // mxFldp = mxGetField(dr, 0,"nspred" );
- const int nsPred = 4; //(int)mxGetScalar(mxFldp);
- // mxFldp = mxGetField(dr, 0,"nboth" );
- const int nBoth = 2; // (int)mxGetScalar(mxFldp);
- // mxFldp = mxGetField(dr, 0,"nfwrd" );
- const int nForw = 5; // 3 (int)mxGetScalar(mxFldp);
- // mxFldp = mxGetField(dr, 0,"nsfwrd" );
- const int nsForw = 7; //(int)mxGetScalar(mxFldp);
-
- // mxFldp = mxGetField(M_, 0,"exo_nbr" );
- const int nExog = 2; // (int)mxGetScalar(mxFldp);
- // mxFldp = mxGetField(M_, 0,"endo_nbr" );
- const int nEndo = 16; // 16(int)mxGetScalar(mxFldp);
- // mxFldp = mxGetField(M_, 0,"param_nbr" );
- const int nPar = 7; //(int)mxGetScalar(mxFldp);
- // it_ should be set to M_.maximum_lag
- // mxFldp = mxGetField(M_, 0,"maximum_lag" );
- const int nMax_lag = 1; //(int)mxGetScalar(mxFldp);
-
- int var_order[] //[18]
- = {
- 5, 6, 8, 10, 11, 12, 14, 7, 13,
- 1, 2, 3, 4, 9, 15, 16
- // 5, 6, 8, 10, 11, 12, 16, 7, 13, 14, 15, 1, 2, 3, 4, 9, 17,
- // 18
- };
- // Vector * varOrder = new Vector(var_order, nEndo);
- vector* var_order_vp = new vector(nEndo); // nEndo));
- for (int v = 0; v < nEndo; v++)
- (*var_order_vp)[v] = var_order[v];
-
- const double ll_incidence[] //[3][18]
- = {1, 5, 21, 2, 6, 22, 0, 7, 23, 0, 8, 24, 0, 9, 0, 0, 10, 0, 3, 11, 0, 0, 12, 0,
- 0, 13, 25, 0, 14, 0, 0, 15, 0, 0, 16, 0, 4, 17, 0, 0, 18, 0, 0, 19, 26, 0, 20, 27};
- TwoDMatrix* llincidence = new TwoDMatrix(3, nEndo, ll_incidence);
-
- const int jcols = nExog + nEndo + nsPred + nsForw; // Num of Jacobian columns
-#ifdef DEBUG
- mexPrintf("k_order_perturbation: jcols= %d .\n", jcols);
-#endif
- // mxFldp= mxGetField(M_, 0,"endo_names" );
- const int nendo = 16; // 16(int)mxGetM(mxFldp);
- const int widthEndo = 6; // (int)mxGetN(mxFldp);
- const char* cNamesCharStr = "mPceWRkdnlggydPc yp A22 __ oo "
- " bb ss ";
- // const char** endoNamesMX= DynareMxArrayToString( mxFldp,nendo,widthEndo);
- const char** endoNamesMX = DynareMxArrayToString(cNamesCharStr, nendo, widthEndo);
-#ifdef DEBUG
- for (int i = 0; i < nEndo; i++)
- {
- mexPrintf("k_ord_perturbation: EndoNameList[%d][0]= %s.\n", i, endoNamesMX[i]);
- }
-#endif
- // mxFldp = mxGetField(M_, 0,"exo_names" );
- const int nexo = 2; //(int)mxGetM(mxFldp);
- const int widthExog = 3; //(int)mxGetN(mxFldp);
- // const char** exoNamesMX= DynareMxArrayToString( mxFldp,nexo,widthExog);
- const char* cExoNamesCharStr = "ee__am";
- const char** exoNamesMX = DynareMxArrayToString(cExoNamesCharStr, nexo, widthExog);
-#ifdef DEBUG
- for (int i = 0; i < nexo; i++)
- {
- mexPrintf("k_ord_perturbation: ExoNameList[%d][0]= %s.\n", i, exoNamesMX[i]);
- }
-#endif
- if ((nEndo != nendo) || (nExog != nexo)) //(nPar != npar)
- {
- mexPrintf("Incorrect number of input parameters.\n");
- return;
- }
-
-#ifdef DEBUG
- for (int i = 0; i < nEndo; i++)
- {
- mexPrintf("k_ord_perturbation: EndoNameList[%d]= %s.\n", i, endoNamesMX[i]);
- }
- // for (int i = 0; i < nPar; i++) {
- //, , mexPrintf("k_ord_perturbation: params_vec[%d]= %g.\n", i, params_vec[i] ); }
- for (int i = 0; i < nPar; i++)
- {
- mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]);
- }
- for (int i = 0; i < nSteady; i++)
- {
- mexPrintf("k_ord_perturbation: ysteady[%d]= %g.\n", i, (*ySteady)[i]);
- }
-
- mexPrintf("k_order_perturbation: nEndo = %d , nExo = %d .\n", nEndo, nExog);
-#endif
- /* Fetch time index */
- // int it_ = (int) mxGetScalar(prhs[3]) - 1;
-
- const int nSteps
- = 0; // Dynare++ solving steps, for time being default to 0 = deterministic steady state
- const double sstol = 1.e-13; // NL solver tolerance from
-
- THREAD_GROUP::max_parallel_threads = 1; // 2 params.num_threads;
-
- try
- {
- // make journal name and journal
- std::string jName(fName); // params.basename);
- jName += ".jnl";
- Journal journal(jName.c_str());
-
-#ifdef DEBUG
- mexPrintf("k_order_perturbation: Call tls init\n");
-#endif
-
- tls.init(kOrder, (nStat + 2 * nPred + 3 * nBoth + 2 * nForw + nExog));
-
-#ifdef DEBUG
- mexPrintf("k_order_perturbation: Calling dynamicDLL constructor.\n");
-#endif
- // DynamicFn * pDynamicFn = loadModelDynamicDLL (fname);
- DynamicModelDLL dynamicDLL(fName, nEndo, jcols, nMax_lag, nExog, dfExt);
-#ifdef DEBUG
- mexPrintf("k_order_perturbation: Calling dynare constructor.\n");
-#endif
- // make KordpDynare object
- KordpDynare dynare(endoNamesMX, nEndo, exoNamesMX, nExog, nPar, // paramNames,
- ySteady, vCov, modParams, nStat, nPred, nForw, nBoth, jcols, NNZD, nSteps,
- kOrder, journal, dynamicDLL, sstol, var_order_vp, // var_order
- llincidence, qz_criterium);
-#ifdef DEBUG
- mexPrintf("k_order_perturbation: Call Approximation constructor \n");
-#endif
- Approximation app(dynare, journal, nSteps, false, qz_criterium);
- // run stochastic steady
-#ifdef DEBUG
- mexPrintf("k_order_perturbation: Calling walkStochSteady.\n");
-#endif
- app.walkStochSteady();
-
- // open mat file
- std::string matfile(fName); //(params.basename);
- matfile += ".mat";
- FILE* matfd = NULL;
- if (NULL == (matfd = fopen(matfile.c_str(), "wb")))
- {
- fprintf(stderr, "Couldn't open %s for writing.\n", matfile.c_str());
- exit(1);
- }
-
-#ifdef DEBUG
- mexPrintf("k_order_perturbation: Filling Mat file outputs.\n");
-#endif
- std::string ss_matrix_name(fName); // params.prefix);
- ss_matrix_name += "_steady_states";
- ConstTwoDMatrix(app.getSS()).writeMat4(matfd, ss_matrix_name.c_str());
-
- // write the folded decision rule to the Mat-4 file
- app.getFoldDecisionRule().writeMat4(matfd, fName); // params.prefix);
-
- fclose(matfd);
-
- map mm;
- app.getFoldDecisionRule().writeMMap(mm, string());
-#ifdef DEBUG
- app.getFoldDecisionRule().print();
- mexPrintf("k_order_perturbation: Map print: \n");
- for (map::const_iterator cit = mm.begin(); cit != mm.end(); ++cit)
- {
- mexPrintf("k_order_perturbation: Map print: string: %s , g:\n", (*cit).first.c_str());
- (*cit).second.print();
- }
-#endif
-
- // get latest ysteady
- double* dYsteady = (dynare.getSteady().base());
- ySteady = (Vector*)(&dynare.getSteady());
- }
- catch (const KordException& e)
- {
- printf("Caugth Kord exception: ");
- e.print();
- return 1; // e.code();
- }
- catch (const TLException& e)
- {
- printf("Caugth TL exception: ");
- e.print();
- return 2; // 255;
- }
- catch (SylvException& e)
- {
- printf("Caught Sylv exception: ");
- e.printMessage();
- return 3; // 255;
- }
- catch (const DynareException& e)
- {
- printf("Caught KordpDynare exception: %s\n", e.message());
- return 4; // 255;
- }
- catch (const ogu::Exception& e)
- {
- printf("Caught ogu::Exception: ");
- e.print();
- return 5; // 255;
- }
-
- // bones for future developement of the Matlab output.
- const int nrhs = 5;
- const int nlhs = 2;
-
- mxArray* prhs[nrhs];
- mxArray* plhs[nlhs];
-
-#ifdef DEBUG
- mexPrintf("k_order_perturbation: Filling MATLAB outputs.\n");
-#endif
-
- double *dgy, *dgu, *ysteady;
- int nb_row_x;
-
- ysteady = NULL;
- if (nlhs >= 1)
- {
- /* Set the output pointer to the output matrix ysteady. */
- plhs[0] = mxCreateDoubleMatrix(nEndo, 1, mxREAL);
- /* Create a C pointer to a copy of the output ysteady. */
- ysteady = mxGetPr(plhs[0]);
- }
-
- dgy = NULL;
- if (nlhs >= 2)
- {
- /* Set the output pointer to the output matrix gy. */
- plhs[1] = mxCreateDoubleMatrix(nEndo, jcols, mxREAL);
- // plhs[1] = (double*)(gy->getData())->base();
- /* Create a C pointer to a copy of the output matrix gy. */
- dgy = mxGetPr(plhs[1]);
- }
-
- dgu = NULL;
- if (nlhs >= 3)
- {
- /* Set the output pointer to the output matrix gu. */
- plhs[2] = mxCreateDoubleMatrix(nEndo, nExog, mxREAL);
- // plhs[2] = (double*)((gu->getData())->base());
- /* Create a C pointer to a copy of the output matrix gu. */
- dgu = mxGetPr(plhs[2]);
- }
-
- return 0;
-}