diff --git a/mex/sources/korderpert/src/k_ord_dynare.cpp b/mex/sources/korderpert/src/k_ord_dynare.cpp
index fcfa7feb8..dc8b07787 100644
--- a/mex/sources/korderpert/src/k_ord_dynare.cpp
+++ b/mex/sources/korderpert/src/k_ord_dynare.cpp
@@ -1,684 +1,675 @@
-/*
- * Copyright (C) 2008-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 <http://www.gnu.org/licenses/>.
- */
-
-// GP, based on work by O.Kamenik
-
-#include <vector>
-#include "first_order.h"
-#include "k_ord_dynare.h"
-
-#include "mex.h" 
-
-#include "memory_file.h"
-
-
-//#include "k_order_perturbation.h"
-#ifndef DYNVERSION
-#define DYNVERSION "unknown"
-#endif
-
-/**************************************************************************************/
-/*       Dynare DynamicModel class                                                                 */
-/**************************************************************************************/
-class KordpJacobian;
-
-KordpDynare::KordpDynare(const char** endo,  int num_endo,
-			   const char** exo, int nexog, int npar, //const char** par,
-   			   Vector* ysteady, TwoDMatrix* vcov, Vector* inParams, int nstat,
-			   int npred, int nforw, int nboth, const int jcols, const int nsteps, int norder, //const char* modName,
-			   Journal& jr, DynamicModelDLL& dynamicDLL, double sstol, const vector<int>* var_order, 
-			   const TwoDMatrix * llincidence, double criterium )
-	: nStat(nstat), nBoth(nboth), nPred(npred), nForw(nforw), nExog(nexog), nPar(npar),
-	nYs(npred + nboth), nYss(nboth + nforw),nY(num_endo), nJcols(jcols), nSteps(nsteps), nOrder(norder), 
-	journal(jr),  dynamicDLL(dynamicDLL), ySteady(ysteady), vCov(vcov), params (inParams), 
-	md(1), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(sstol), varOrder( var_order), 
-	ll_Incidence(llincidence), qz_criterium(criterium) 
-{
-#ifdef DEBUG		
-   mexPrintf("k_ord_dynare Dynare constructor: ny=%d, order=%d, nPar=%d .\n", nY,nOrder,nPar);
-	for (int i = 0; i < nY; i++) {
-        mexPrintf("k_ord_dynare calling DynareNameList names[%d]= %s.\n", i, endo[i] );}
-	for (int i = 0; i < nPar; i++) {
-        mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*params)[i]);  }
-	for (int i = 0; i < nY; i++) {
-        mexPrintf("k_ord_perturbation: ysteady[%d]= %g.\n", i, (*ySteady)[i]);  }
-    mexPrintf("k_ord_dynare: dynare constructor, trying namelists.\n");
-#endif		
-	try{
-		dnl = new DynareNameList(*this, endo);
-		denl = new DynareExogNameList(*this, exo);
-#ifdef DEBUG		
-    mexPrintf("k_ord_dynare: dynare constructor, trying StateNamelist.\n");
-#endif		
-		dsnl = new DynareStateNameList(*this, *dnl, *denl);
-
-		JacobianIndices= ReorderDynareJacobianIndices( varOrder);
-
-	//	Initialise ModelDerivativeContainer(*this, this->md, nOrder);
-    	for (int iord = 1; iord <= nOrder; iord++) {
-    		FSSparseTensor* t = new FSSparseTensor(iord, nY+nYs+nYss+nExog,nY);
-        	md.insert(t);
-        }
-
-	}
-	catch (...){
-        mexPrintf("k_ord_dynare: dynare constructor, error in StateNamelist construction.\n");
-        throw DynareException(__FILE__, __LINE__, string("Could not construct Name Lists. \n"));
-    }
-}
-
-KordpDynare::KordpDynare(const KordpDynare& dynare)
-	: nStat(dynare.nStat), nBoth(dynare.nBoth),	nPred(dynare.nPred), 
-	nForw(dynare.nForw), nExog(dynare.nExog),  nPar(dynare.nPar),
-	nYs(dynare.nYs), nYss(dynare.nYss),nY(dynare.nY), nJcols(dynare.nJcols), 
-	nSteps(dynare.nSteps), nOrder(dynare.nOrder), journal(dynare.journal),
-	dynamicDLL(dynare.dynamicDLL), //modName(dynare.modName),
-	ySteady(NULL), params(NULL), vCov(NULL), md(dynare.md), 
-	dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(dynare.ss_tol), 
-	varOrder(dynare.varOrder), ll_Incidence(dynare.ll_Incidence),
-	JacobianIndices(dynare.JacobianIndices), qz_criterium(dynare.qz_criterium)
-{
-	ySteady = new Vector(*(dynare.ySteady));
-	params = new Vector(*(dynare.params));
-	vCov = new TwoDMatrix(*(dynare.vCov));
-	dnl = new DynareNameList(dynare);//(*this);
-	denl = new DynareExogNameList(dynare);//(*this);
-	dsnl = new DynareStateNameList(*this, *dnl, *denl);
-}
-
-KordpDynare::~KordpDynare()
-{
-	if (ySteady)
-		delete ySteady;
-	if (params)
-		delete params;
-	if (vCov)
-		delete vCov;
-	if (dnl)
-		delete dnl;
-	if (dsnl)
-		delete dsnl;
-	if (denl)
-		delete denl;
-}
-
-
-/** This clears the container of model derivatives and initializes it
- * inserting empty sparse tensors up to the given order. */
-ModelDerivativeContainer::ModelDerivativeContainer(const KordpDynare& model, 
-		TensorContainer<FSSparseTensor>& mod_ders, int order): md(mod_ders)
-{
-	md.clear();
-	for (int iord = 1; iord <= order; iord++) {
-		FSSparseTensor* t = new FSSparseTensor(iord, model.ny()+model.nys()+model.nyss()+model.nexog(), model.ny());
-		md.insert(t);
-	}
-}
-
-
-void KordpDynare::solveDeterministicSteady(Vector& steady)
-{
-	JournalRecordPair pa(journal);
-	pa << "Non-linear solver for deterministic steady state By-passed " << endrec;
-	/*************************;  GP Dec 08 by-pass
-	KordpVectorFunction dvf(*this);
-	KordpJacobian dj(*this);
-	ogu::NLSolver nls(dvf, dj, 500, ss_tol, journal);
-	int iter;
-	if (! nls.solve(*ySteady, iter))
-		throw DynareException(__FILE__, __LINE__,
-							  "Could not obtain convergence in non-linear solver");
-	***************************/
-}
-
-// evaluate system at given y_t=y_{t+1}=y_{t-1}, and given shocks x_t
-void KordpDynare::evaluateSystem(Vector& out, const Vector& yy, const Vector& xx)
-{
-	dynamicDLL.eval( yy,  xx, //int nb_row_x, 
-				 params, //int it_, 
-                 out, NULL, NULL);
-
-}
-
-// evaluate system at given y^*_{t-1}, y_t, y^{**}_{t+1} and at
-// exogenous x_t, all three vectors yym, yy, and yyp have the
-// respective lengths of y^*_{t-1}, y_t, y^{**}_{t+1}
-
-void KordpDynare::evaluateSystem(Vector& out, const Vector& yym, const Vector& yy,
-							const Vector& yyp, const Vector& xx)
-{
-#ifdef DEBUG
-	mexPrintf("k_order_dynaare.cpp: Call eval in EvaluateSystem\n");
-#endif
-	dynamicDLL.eval( yy,  xx, //int nb_row_x, 
-				 params, //int it_, 
-                 out, NULL, NULL);
-}
-/************************************************
-* this is main derivative calculation functin that indirectly calls dynamic.dll
-* which performs actual calculation and reorders
-***************************************************/ 
-//void KordpDynare::calcDerivatives(const Vector& yy, const TwoDMatrix& xx)
-void KordpDynare::calcDerivatives(const Vector& yy, const Vector& xx)
-{
-   // Hessian TwoDMatrix *g2;
-	TwoDMatrix *g2=NULL;
-    //Jacobian 
-	TwoDMatrix * g1=new TwoDMatrix(nY,nJcols); // generate g1 for jacobian  
-	g1->zeros();
-
-	if ((nJcols != g1->ncols()) && ( nY != g1->nrows())) {
-		mexPrintf("k_ord_dynare.cpp: Error in calcDerivatives: Created wrong jacobian");
-		return;
-	}
-
-    // Hessian TwoDMatrix *g2;
-	if (nOrder>1){
-	//TwoDMatrix * 
-		g2=new TwoDMatrix(nY,nJcols*nJcols); // generate g2 for Hessian  
-		g2->zeros();
-    }
-	Vector& out= *(new Vector(nY));
-	out.zeros();
-	const Vector * llxYYp; // getting around the constantness
-	if ((nJcols - nExog) > yy.length()){
-		llxYYp=  (LLxSteady( yy));
-	} else {
-		llxYYp= &yy;
-	}
-	const Vector & llxYY=*(llxYYp);
-
-#ifdef DEBUG
-	mexPrintf("k_order_dynaare.cpp: Call eval in calcDerivatives\n");
-#endif
-	try {
-			dynamicDLL.eval( llxYY,  xx, //int nb_row_x, 
-					params, //int it_, 
-					out, g1, g2);
-//		}
-	}
-	catch (...){
-			mexPrintf("k_ord_dynare.cpp: Error in dynamicDLL.eval in calcDerivatives");
-			return;
-	}
-	if ((nJcols!=g1->ncols()) && ( nY != g1->nrows())) {
-			mexPrintf("k_ord_dynare.cpp: Error in calcDerivatives: dynamicDLL.eval returned wrong jacobian");
-			return;
-	}
-	//	ReorderCols(g1, JacobianIndices); and populate container
-	populateDerivativesContainer(g1,1,JacobianIndices);
-	if (nOrder>1){
-	  //		ReorderBlocks(g2,JacobianIndices);
-	  populateDerivativesContainer(g2,2,JacobianIndices);
-	}
-
-}
-/* This version is not currently in use */
-void KordpDynare::calcDerivatives(const Vector& yy, ogu::Jacobian& jacob)
-{
-//	ConstVector yym(yy, nstat(), nys());
-//	ConstVector yyp(yy, nstat()+npred(), nyss());
-
-//	Vector yyp(yy, nstat()+npred(), nyss());
-
-	//double *g1, *g2;
-	TwoDMatrix * jj= &jacob;
-    Vector& out= *(new Vector(nY)); 
-	out.zeros();
-    Vector& xx= *(new Vector(nExog)); 
-	xx.zeros();
-	dynamicDLL.eval( yy,  xx, //int nb_row_x, 
-				params, //int it_, 
-				out, jj, NULL);
-   //    model derivatives FSSparseTensor instance
-        FSSparseTensor &mdTi=*(new FSSparseTensor (1, jj->ncols(),jj->nrows())); 
-        for (int i = 0; i<jj->ncols(); i++){
-                for (int j = 0; j<jj->nrows(); j++){
-                    if (jj->get(i,j)!=0.0) // populate sparse if not zero
-                        mdTi.insert(i, j, jj->get(i,j));
-                }
-        }
-        // md container
-        md.insert(&mdTi);
-		delete &out;
-		delete &xx;
-}
-
-void KordpDynare::calcDerivativesAtSteady()
-{
-	Vector xx(nexog());
-	xx.zeros();
-	calcDerivatives(*ySteady, xx);
-}
-
-/*******************************************************************************
-* populateDerivatives to sparse Tensor and fit it in the Derivatives Container
-*******************************************************************************/
-void KordpDynare::populateDerivativesContainer(TwoDMatrix*g, int ord, const vector<int>*  vOrder)
-{
-
-#ifdef DEBUG
-	mexPrintf("k_ord_dynare.cpp: populate FSSparseTensor in calcDerivatives: cols=%d , rows=%d\n"
-        , g->ncols(),g->nrows());
-#endif
-
-	// model derivatives FSSparseTensor instance 
-    FSSparseTensor *mdTi=(new FSSparseTensor (ord, nJcols,g->nrows()));
-
-    IntSequence s(ord,0);
-    s[0] = 0;
-    s[1] = 0;
-    int i = 0;
-    while (i < g->ncols()){
-
-      // insert new elements in each row
-      if (ord == 1){
-		for (int j = 0; j < g->nrows(); j++){
-		  double x;
-		  if (s[0] < nJcols-nExog)
-			x = g->get(j,(*vOrder)[s[0]]);
-		  else
-			x = g->get(j,s[0]);
-		  if (x != 0.0)
-			mdTi->insert(s, j, x);
-		}
-		s[0]++;
-      }
-      else{	
-		int s0, s1;
-		if (s[0] < nJcols-nExog)
-		  s0 = (*vOrder)[s[0]];
-		else
-		  s0 = s[0];
-		if (s[1] < nJcols-nExog)
-		  s1 = (*vOrder)[s[1]];
-		else
-		  s1 = s[1];
-		if (s[1] >= s[0]){
-		  s.print();
-		  std::cout << s0 << " " << s1 << "\n";
-		  int i1 = s0*nJcols+s1;
-		  for (int j = 0; j < g->nrows(); j++){
-			double x = g->get(j,i1);
-			if (x != 0.0)
-			  mdTi->insert(s, j, x);
-		  }
-		}
-		s[1]++;
-		// when one order index is finished
-		// increase the previous one
-		if (s[1] == nJcols){
-		  s[0]++;
-		  // update starting position of next indices
-		  // in order to avoid dealing twice with the same
-		  // symmetry. Increase matrix column counter
-		  // accordingly
-		  s[1] = 0;
-		}
-      }
-
-      i++;
-
-    }
-    mdTi->print();
-    // md container
-    //md.clear();// this is to be used only for 1st order!!
-    md.remove(Symmetry(ord));
-    md.insert(mdTi);//(&mdTi);
-}
-
-
-void KordpDynare::writeModelInfo(Journal& jr) const
-{
-	// write info on variables
-	{
-		JournalRecordPair rp(journal);
-		rp << "Information on variables" << endrec;
-		JournalRecord rec1(journal);
-		rec1 << "Number of endogenous:            " << ny() << endrec;
-		JournalRecord rec2(journal);
-		rec2 << "Number of exogenous:             " << nexog() << endrec;
-		JournalRecord rec3(journal);
-		rec3 << "Number of static:                " << nstat() << endrec;
-		JournalRecord rec4(journal);
-		rec4 << "Number of predetermined:         " << npred()+nboth() << endrec;
-		JournalRecord rec5(journal);
-		rec5 << "Number of forward looking:       " << nforw()+nboth() << endrec;
-		JournalRecord rec6(journal);
-		rec6 << "Number of both:                  " << nboth() << endrec;
-	}
-
-}
-/*********************************************************
-* LLxSteady()
-* returns ySteady extended with leads and lags suitable for 
-* passing to <model>_dynamic DLL
-*************************************************************/
-Vector * KordpDynare::LLxSteady( const Vector& yS){
-	if ((nJcols-nExog) == yS.length()) {
-		mexPrintf("k_ord_dynare.cpp: Warning in LLxSteady: ySteady already. right size");
-		return NULL;
-	}
-	// create temporary square 2D matrix size nEndo x nEndo (sparse)  
-	// for the lag, current and lead blocks of the jacobian
-	Vector * llxSteady = new Vector(nJcols-nExog); 
-	try{
-		for (int ll_row=0; ll_row< ll_Incidence->nrows(); ll_row++)
-		{
-			// populate (non-sparse) vector with ysteady values 
-			for (int i=0;i<nY;i++){
-				if(ll_Incidence->get(ll_row,i))
-					(*llxSteady)[((int)ll_Incidence->get(ll_row,i))-1] = yS[i];
-			}
-		}
-	} catch (const TLException& e) {
-		mexPrintf("Caugth TL exception in LLxSteady: ");
-		e.print();
-		return NULL;// 255;
-	}catch (...){
-		mexPrintf(" Error in LLxSteady - wrong index?");
-	}
-
-#ifdef DEBUG		
-	for (int j=0;j<nJcols-nExog;j++)
-			mexPrintf("LLxSteady: [%d] =%f .\n", j, (*llxSteady)[j]);
-#endif	
-	return llxSteady;
-}
-
-
-/************************************
-* Reorder DynareJacobianIndices of variables in a vector according to
-* given int * varOrder together with lead & lag incidence matrix and 
-* any the extra columns for exogenous vars, and then, 
-* reorders its blocks given by the varOrder and the Dynare++ expectations:
-
-* extra	nboth+ npred (t-1) lags
-* varOrder
-		static:
-		pred
-		both    
-		forward    
-* extra both + nforw (t+1) leads, and 
-* extra exogen
-
-* so to match the jacobian organisation expected by the Appoximation class 
-		both + nforw (t+1) leads
-		static
-		pred
-		both
-		forward
-		nboth+ npred  (t-1) lags
-		exogen
-************************************/
-
-vector<int> * KordpDynare::ReorderDynareJacobianIndices( const vector<int> * varOrder){
-//	if ((nJcols != tdx->ncols()) && ( nY != tdx->nrows())) {
-//		mexPrintf("k_ord_dynare.cpp: Error in ReorderBlocks: wrong size of jacobian");
-//		return;
-//	}
-	// create temporary square 2D matrix size nEndo x nEndo (sparse)  
-	// for the lag, current and lead blocks of the jacobian
-//	int * JacobianIndices = (int*) calloc(nJcols+1, sizeof(int)); 
-	vector<int> * JacobianIndices = new vector<int>(nJcols); 
-	vector <int> tmp(nY); 
-	int i,j, rjoff=nJcols-nExog-1; //, ll_off, j;
-#ifdef DEBUG		
-	mexPrintf("ReorderDynareJacobianIndice:ll_Incidence->nrows() =%d .\n", ll_Incidence->nrows());
-#endif
-	try{
-		for (int ll_row=0; ll_row< ll_Incidence->nrows(); ll_row++)
-		{
-			// reorder in orde-var order & populate temporary nEndo (sparse) vector with 
-			// the lag, current and lead blocks of the jacobian respectively
-			for (i=0;i<nY;i++){
-//				tmp[varOrder[j]]=(int)ll_Incidence->get(ll_row,j); 
-				tmp[i]=((int)ll_Incidence->get(ll_row,(*varOrder)[i]-1)); 
-#ifdef DEBUG		
-	mexPrintf("get(ll_row,(*varOrder)[%d]-1)) = tmp[%d]=%d .\n", 
-		i, i, (int)ll_Incidence->get(ll_row,(*varOrder)[i]-1));
-#endif
-			}
-			// write the reordered blocks back to the jacobian
-			// in reverse order
-			for (j=nY-1;j>=0;j--){
-				if (tmp[j]){
-					(*JacobianIndices)[rjoff]=tmp[j] -1;
-					rjoff--;
-					if (rjoff<0){
-	//					mexPrintf(" Error in ReorderIndices - negative rjoff index?");
-						break;
-	//					return NULL;
-					}
-				}
-			}
-		}
-	} catch (const TLException& e) {
-		mexPrintf("Caugth TL exception in ReorderIndices: ");
-		e.print();
-		return NULL;// 255;
-	}catch (...){
-		mexPrintf(" Error in ReorderIndices - wrong index?");
-	}
-	//add the indices for the nExog exogenous jacobians
-	for (j=nJcols-nExog;j<nJcols;j++){
-		(*JacobianIndices)[j]=j;
-	}
-#ifdef DEBUG		
-	for (j=0;j<nJcols;j++)
-			mexPrintf("ReorderDynareJacobianIndice: [%d] =%d .\n", j, (*JacobianIndices)[j]);
-#endif	
-	return JacobianIndices;
-}
-
-
-/************************************
-* Reorder first set of columns of variables in a (jacobian) matrix 
-* according to order given in  varsOrder together with the extras 
-* assuming tdx ncols() - nExog is eaqual or less than length of varOrder and
-* of any of its elements too.
-************************************/
-
-void KordpDynare::ReorderCols(TwoDMatrix * tdx, const vector<int> * vOrder){
-
-	if (tdx->ncols() > vOrder->size()){
-		mexPrintf(" Error in ReorderColumns - size of order var is too small");
-		return;
-	}
-	TwoDMatrix tmp(*tdx); // temporary 2D matrix
-	TwoDMatrix &tmpR=tmp;
-	tdx->zeros();// empty original matrix
-	// reorder the columns
-	try{
-		for (int i =0; i<tdx->ncols() ; i++)
-			tdx->copyColumn(tmpR,(*vOrder)[i],i);
-	} catch (const TLException& e) {
-		printf("Caugth TL exception in ReorderColumns: ");
-		e.print();
-		return;// 255;
-	}catch (...){
-		mexPrintf(" Error in ReorderColumns - wrong index?");
-	}
-}
-void KordpDynare::ReorderCols(TwoDMatrix * tdx, const int * vOrder){
-
-	TwoDMatrix tmp(*tdx); // temporary 2D matrix
-	TwoDMatrix &tmpR=tmp;
-	tdx->zeros();// empty original matrix
-	// reorder the columns
-	try{
-		for (int i =0; i<tdx->ncols() ; i++)
-			tdx->copyColumn(tmpR,vOrder[i],i);
-	} catch (const TLException& e) {
-		printf("Caugth TL exception in ReorderColumns: ");
-		e.print();
-		return;// 255;
-	}catch (...){
-		mexPrintf(" Error in ReorderColumns - wrong index?");
-	}
-}
-
-/***********************************************************************
-* Recursive hierarchical block reordering of the higher order, input model 
-*	derivatives inc. Hessian   
-* This is now obsolete but kept in in case it is needed
-***********************************************************************/
-
-void KordpDynare::ReorderBlocks(TwoDMatrix * tdx, const vector<int> * vOrder){
-	// determine order of the matrix
-
-	double dbOrder = log(tdx->ncols())/log(nJcols);
-	int ibOrder= (int) dbOrder;
-	if ((double )ibOrder != dbOrder || ibOrder>nOrder) {
-		mexPrintf(" Error in ReorderBlocks - wrong order %d", dbOrder);
-		return;
-	}
-	TwoDMatrix tmp(*tdx); // temporary 2D matrix
-	TwoDMatrix &tmpR=tmp;
-	tdx->zeros();// empty original matrix
-
-	if(ibOrder>1){
-		int nBlocks=tmp.ncols()/ nJcols; //pow((float)nJcols,ibOrder-1);
-		int bSize=tmp.ncols()/nBlocks;
-		for (int j = 0; j<nBlocks;  ++j){
-			TwoDMatrix subtdx(tmpR, bSize*((*vOrder)[j]), bSize);
-			ReorderBlocks(&subtdx, vOrder);
-			tdx->place(subtdx, 0, bSize*j);
-		}
-	} else{
-		//ReorderColumns(TwoDMatrix * subtdx, const vector<int> * vOrder)
-		if (tdx->ncols() > vOrder->size()){
-			mexPrintf(" Error in ReorderColumns - size of order var is too small");
-			return;
-		}
-		// reorder the columns
-		try{
-			for (int i =0; i<tdx->ncols() ; i++)
-				tdx->copyColumn(tmpR,(*vOrder)[i],i);
-		} catch (const TLException& e) {
-			printf("Caugth TL exception in ReorderColumns: ");
-			e.print();
-			return;// 255;
-		}catch (...){
-			mexPrintf(" Error in ReorderColumns - wrong index?");
-		}
-	}
-}
-
-
-/****************************
-*  K-Order Perturbation instance of Jacobian:
-************************************/
-KordpJacobian::KordpJacobian(KordpDynare& dyn)
-	: Jacobian(dyn.ny()), dyn(dyn)
-{
-	zeros();
-};
-
-void KordpJacobian::eval(const Vector& yy)
-{
-		dyn.calcDerivatives( yy, *this);
-
-};
-
-void KordpVectorFunction::eval(const ConstVector& in, Vector& out)
-{
-	check_for_eval(in, out);
-	Vector xx(d.nexog());
-	xx.zeros();
-	d.evaluateSystem(out, in, xx);
-}
-
-/**************************************************************************************/
-/*       DynareNameList class                                                         */
-/**************************************************************************************/
-vector<int> DynareNameList::selectIndices(const vector<const char*>& ns) const
-{
-	vector<int> res;
-	for (unsigned int i = 0; i < ns.size(); i++) {
-		int j = 0;
-		while (j < getNum() && strcmp(getName(j), ns[i]) != 0)
-			j++;
-		if (j == getNum())
-			throw DynareException(__FILE__, __LINE__,
-								  string("Couldn't find name for ") + ns[i] +
-								  " in DynareNameList::selectIndices");
-		res.push_back(j);
-	}
-	return res;
-}
-
-DynareNameList::DynareNameList(const  KordpDynare& dynare)
-{
-	for (int i = 0; i < dynare.ny(); i++) {
-		names.push_back(dynare.dnl->getName(i));
-	}
-}
-DynareNameList::DynareNameList(const KordpDynare& dynare, const char ** namesp)
-{
-#ifdef DEBUG		
-    mexPrintf("k_ord_dynare DynareNameList.\n");
-    mexPrintf("k_ord_dynare DynareNameList dynare.ny=%d .\n", dynare.ny());
-#endif		
-    
-	for (int i = 0; i < dynare.ny(); i++) {
-#ifdef DEBUG		
-    mexPrintf("k_ord_dynare DynareNameList names[%d]= %s.\n", i, namesp[i] );
-#endif		
-		names.push_back(namesp[i]);
-	}
-}
-
-DynareExogNameList::DynareExogNameList(const KordpDynare& dynare)
-{
-	for (int i = 0; i < dynare.nexog(); i++) {
-		names.push_back(dynare.denl->getName(i));
-	}
-}
-
-DynareExogNameList::DynareExogNameList(const KordpDynare& dynare, const char ** namesp)
-{
-#ifdef DEBUG		
-    mexPrintf("k_ord_dynare DynareExogNameList dynare.nexog=%d .\n", dynare.nexog());
-#endif		
-	for (int i = 0; i < dynare.nexog(); i++) {
-#ifdef DEBUG		
-    mexPrintf("k_ord_dynare DynareExogNameList names[%d]= %s.\n", i, namesp[i] );
-#endif		
-		names.push_back(namesp[i]);
-	}
-}
-
-DynareStateNameList::DynareStateNameList(const KordpDynare& dynare, const DynareNameList& dnl,
-										 const DynareExogNameList& denl)
-{
-	for (int i = 0; i < dynare.nys(); i++){
-#ifdef DEBUG		
-    mexPrintf("k_ord_dynare DynareStateNameList dnl names[%d]= %s.\n", i, dnl.getName(i+dynare.nstat()) );
-#endif		
-		names.push_back(dnl.getName(i+dynare.nstat()));
-    }
-	for (int i = 0; i < dynare.nexog(); i++){
-#ifdef DEBUG		
-    mexPrintf("k_ord_dynare DynareStateNameList denl names[%d]= %s.\n", i, denl.getName(i));
-#endif		
-		names.push_back(denl.getName(i));
-    }
-}
-
+/*
+* Copyright (C) 2008-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 <http://www.gnu.org/licenses/>.
+*/
+
+// GP, based on work by O.Kamenik
+
+#include <vector>
+#include "first_order.h"
+#include "k_ord_dynare.h"
+
+#include "mex.h" 
+
+#include "memory_file.h"
+
+
+//#include "k_order_perturbation.h"
+#ifndef DYNVERSION
+#define DYNVERSION "unknown"
+#endif
+
+/**************************************************************************************/
+/*       Dynare DynamicModel class                                                                 */
+/**************************************************************************************/
+class KordpJacobian;
+
+KordpDynare::KordpDynare(const char** endo,  int num_endo,
+             const char** exo, int nexog, int npar, //const char** par,
+             Vector* ysteady, TwoDMatrix* vcov, Vector* inParams, int nstat,
+             int npred, int nforw, int nboth, const int jcols, const int nsteps, int norder, //const char* modName,
+             Journal& jr, DynamicModelDLL& dynamicDLL, double sstol, const vector<int>* var_order, 
+             const TwoDMatrix * llincidence, double criterium )
+             : nStat(nstat), nBoth(nboth), nPred(npred), nForw(nforw), nExog(nexog), nPar(npar),
+             nYs(npred + nboth), nYss(nboth + nforw),nY(num_endo), nJcols(jcols), nSteps(nsteps), nOrder(norder), 
+             journal(jr),  dynamicDLL(dynamicDLL), ySteady(ysteady), vCov(vcov), params (inParams), 
+             md(1), dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(sstol), varOrder( var_order), 
+             ll_Incidence(llincidence), qz_criterium(criterium) 
+{
+#ifdef DEBUG		
+  mexPrintf("k_ord_dynare Dynare constructor: ny=%d, order=%d, nPar=%d .\n", nY,nOrder,nPar);
+  for (int i = 0; i < nY; i++) {
+    mexPrintf("k_ord_dynare calling DynareNameList names[%d]= %s.\n", i, endo[i] );}
+  for (int i = 0; i < nPar; i++) {
+    mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*params)[i]);  }
+  for (int i = 0; i < nY; i++) {
+    mexPrintf("k_ord_perturbation: ysteady[%d]= %g.\n", i, (*ySteady)[i]);  }
+  mexPrintf("k_ord_dynare: dynare constructor, trying namelists.\n");
+#endif		
+  try{
+    dnl = new DynareNameList(*this, endo);
+    denl = new DynareExogNameList(*this, exo);
+#ifdef DEBUG		
+    mexPrintf("k_ord_dynare: dynare constructor, trying StateNamelist.\n");
+#endif		
+    dsnl = new DynareStateNameList(*this, *dnl, *denl);
+    
+    JacobianIndices= ReorderDynareJacobianIndices( varOrder);
+    
+    //	Initialise ModelDerivativeContainer(*this, this->md, nOrder);
+    for (int iord = 1; iord <= nOrder; iord++) {
+    		FSSparseTensor* t = new FSSparseTensor(iord, nY+nYs+nYss+nExog,nY);
+        md.insert(t);
+    }
+  }
+  catch (...){
+    mexPrintf("k_ord_dynare: dynare constructor, error in StateNamelist construction.\n");
+    throw DynareException(__FILE__, __LINE__, string("Could not construct Name Lists. \n"));
+  }
+}
+
+KordpDynare::KordpDynare(const KordpDynare& dynare)
+      : nStat(dynare.nStat), nBoth(dynare.nBoth),	nPred(dynare.nPred), 
+      nForw(dynare.nForw), nExog(dynare.nExog),  nPar(dynare.nPar),
+      nYs(dynare.nYs), nYss(dynare.nYss),nY(dynare.nY), nJcols(dynare.nJcols), 
+      nSteps(dynare.nSteps), nOrder(dynare.nOrder), journal(dynare.journal),
+      dynamicDLL(dynare.dynamicDLL), //modName(dynare.modName),
+      ySteady(NULL), params(NULL), vCov(NULL), md(dynare.md), 
+      dnl(NULL), denl(NULL), dsnl(NULL), ss_tol(dynare.ss_tol), 
+      varOrder(dynare.varOrder), ll_Incidence(dynare.ll_Incidence),
+      JacobianIndices(dynare.JacobianIndices), qz_criterium(dynare.qz_criterium)
+{
+  ySteady = new Vector(*(dynare.ySteady));
+  params = new Vector(*(dynare.params));
+  vCov = new TwoDMatrix(*(dynare.vCov));
+  dnl = new DynareNameList(dynare);//(*this);
+  denl = new DynareExogNameList(dynare);//(*this);
+  dsnl = new DynareStateNameList(*this, *dnl, *denl);
+}
+
+KordpDynare::~KordpDynare()
+{
+  if (ySteady)
+    delete ySteady;
+  if (params)
+    delete params;
+  if (vCov)
+    delete vCov;
+  if (dnl)
+    delete dnl;
+  if (dsnl)
+    delete dsnl;
+  if (denl)
+    delete denl;
+}
+
+
+/** This clears the container of model derivatives and initializes it
+* inserting empty sparse tensors up to the given order. */
+ModelDerivativeContainer::ModelDerivativeContainer(const KordpDynare& model, 
+                                                   TensorContainer<FSSparseTensor>& mod_ders, int order): md(mod_ders)
+{
+  md.clear();
+  for (int iord = 1; iord <= order; iord++) {
+    FSSparseTensor* t = new FSSparseTensor(iord, model.ny()+model.nys()+model.nyss()+model.nexog(), model.ny());
+    md.insert(t);
+  }
+}
+
+
+void KordpDynare::solveDeterministicSteady(Vector& steady)
+{
+  JournalRecordPair pa(journal);
+  pa << "Non-linear solver for deterministic steady state By-passed " << endrec;
+  /*************************;  GP Dec 08 by-pass
+  KordpVectorFunction dvf(*this);
+  KordpJacobian dj(*this);
+  ogu::NLSolver nls(dvf, dj, 500, ss_tol, journal);
+  int iter;
+  if (! nls.solve(*ySteady, iter))
+		throw DynareException(__FILE__, __LINE__,
+    "Could not obtain convergence in non-linear solver");
+  ***************************/
+}
+
+// evaluate system at given y_t=y_{t+1}=y_{t-1}, and given shocks x_t
+void KordpDynare::evaluateSystem(Vector& out, const Vector& yy, const Vector& xx)
+{
+  dynamicDLL.eval( yy,  xx, //int nb_row_x, 
+    params, //int it_, 
+    out, NULL, NULL);
+  
+}
+
+// evaluate system at given y^*_{t-1}, y_t, y^{**}_{t+1} and at
+// exogenous x_t, all three vectors yym, yy, and yyp have the
+// respective lengths of y^*_{t-1}, y_t, y^{**}_{t+1}
+
+void KordpDynare::evaluateSystem(Vector& out, const Vector& yym, const Vector& yy,
+                                 const Vector& yyp, const Vector& xx)
+{
+#ifdef DEBUG
+  mexPrintf("k_order_dynaare.cpp: Call eval in EvaluateSystem\n");
+#endif
+  dynamicDLL.eval( yy,  xx, //int nb_row_x, 
+    params, //int it_, 
+    out, NULL, NULL);
+}
+/************************************************
+* this is main derivative calculation functin that indirectly calls dynamic.dll
+* which performs actual calculation and reorders
+***************************************************/ 
+//void KordpDynare::calcDerivatives(const Vector& yy, const TwoDMatrix& xx)
+void KordpDynare::calcDerivatives(const Vector& yy, const Vector& xx)
+{
+  // Hessian TwoDMatrix *g2;
+  TwoDMatrix *g2=NULL;
+  //Jacobian 
+  TwoDMatrix * g1=new TwoDMatrix(nY,nJcols); // generate g1 for jacobian  
+  g1->zeros();
+  
+  if ((nJcols != g1->ncols()) && ( nY != g1->nrows())) {
+    mexPrintf("k_ord_dynare.cpp: Error in calcDerivatives: Created wrong jacobian");
+    return;
+  }
+  
+  // Hessian TwoDMatrix *g2;
+  if (nOrder>1){
+    //TwoDMatrix * 
+    g2=new TwoDMatrix(nY,nJcols*nJcols); // generate g2 for Hessian  
+    g2->zeros();
+  }
+  Vector& out= *(new Vector(nY));
+  out.zeros();
+  const Vector * llxYYp; // getting around the constantness
+  if ((nJcols - nExog) > yy.length()){
+    llxYYp=  (LLxSteady( yy));
+  } else {
+    llxYYp= &yy;
+  }
+  const Vector & llxYY=*(llxYYp);
+  
+#ifdef DEBUG
+  mexPrintf("k_order_dynaare.cpp: Call eval in calcDerivatives\n");
+#endif
+  try {
+    dynamicDLL.eval( llxYY,  xx, //int nb_row_x, 
+      params, //int it_, 
+      out, g1, g2);
+    //		}
+  }
+  catch (...){
+    mexPrintf("k_ord_dynare.cpp: Error in dynamicDLL.eval in calcDerivatives");
+    return;
+  }
+  if ((nJcols!=g1->ncols()) && ( nY != g1->nrows())) {
+    mexPrintf("k_ord_dynare.cpp: Error in calcDerivatives: dynamicDLL.eval returned wrong jacobian");
+    return;
+  }
+  //	ReorderCols(g1, JacobianIndices); and populate container
+  populateDerivativesContainer(g1,1,JacobianIndices);
+  if (nOrder>1){
+    //		ReorderBlocks(g2,JacobianIndices);
+    populateDerivativesContainer(g2,2,JacobianIndices);
+  }
+  
+}
+/* This version is not currently in use */
+void KordpDynare::calcDerivatives(const Vector& yy, ogu::Jacobian& jacob)
+{
+  
+  //double *g1, *g2;
+  TwoDMatrix * jj= &jacob;
+  Vector& out= *(new Vector(nY)); 
+  out.zeros();
+  Vector& xx= *(new Vector(nExog)); 
+  xx.zeros();
+  dynamicDLL.eval( yy,  xx, //int nb_row_x, 
+				params, //int it_, 
+        out, jj, NULL);
+  //    model derivatives FSSparseTensor instance
+  FSSparseTensor &mdTi=*(new FSSparseTensor (1, jj->ncols(),jj->nrows())); 
+  for (int i = 0; i<jj->ncols(); i++){
+    for (int j = 0; j<jj->nrows(); j++){
+      if (jj->get(i,j)!=0.0) // populate sparse if not zero
+        mdTi.insert(i, j, jj->get(i,j));
+    }
+  }
+  // md container
+  md.insert(&mdTi);
+		delete &out;
+    delete &xx;
+}
+
+void KordpDynare::calcDerivativesAtSteady()
+{
+  Vector xx(nexog());
+  xx.zeros();
+  calcDerivatives(*ySteady, xx);
+}
+
+/*******************************************************************************
+* populateDerivatives to sparse Tensor and fit it in the Derivatives Container
+*******************************************************************************/
+void KordpDynare::populateDerivativesContainer(TwoDMatrix*g, int ord, const vector<int>*  vOrder)
+{
+  
+#ifdef DEBUG
+  mexPrintf("k_ord_dynare.cpp: populate FSSparseTensor in calcDerivatives: cols=%d , rows=%d\n"
+    , g->ncols(),g->nrows());
+#endif
+  
+  // model derivatives FSSparseTensor instance 
+  FSSparseTensor *mdTi=(new FSSparseTensor (ord, nJcols,g->nrows()));
+  
+  IntSequence s(ord,0);
+  s[0] = 0;
+  s[1] = 0;
+  int i = 0;
+  while (i < g->ncols()){
+    
+    // insert new elements in each row
+    if (ord == 1){
+      for (int j = 0; j < g->nrows(); j++){
+        double x;
+        if (s[0] < nJcols-nExog)
+          x = g->get(j,(*vOrder)[s[0]]);
+        else
+          x = g->get(j,s[0]);
+        if (x != 0.0)
+          mdTi->insert(s, j, x);
+      }
+      s[0]++;
+    }
+    else{	
+      int s0, s1;
+      if (s[0] < nJcols-nExog)
+        s0 = (*vOrder)[s[0]];
+      else
+        s0 = s[0];
+      if (s[1] < nJcols-nExog)
+        s1 = (*vOrder)[s[1]];
+      else
+        s1 = s[1];
+      if (s[1] >= s[0]){
+        s.print();
+        std::cout << s0 << " " << s1 << "\n";
+        int i1 = s0*nJcols+s1;
+        for (int j = 0; j < g->nrows(); j++){
+          double x = g->get(j,i1);
+          if (x != 0.0)
+            mdTi->insert(s, j, x);
+        }
+      }
+      s[1]++;
+      // when one order index is finished
+      // increase the previous one
+      if (s[1] == nJcols){
+        s[0]++;
+        // update starting position of next indices
+        // in order to avoid dealing twice with the same
+        // symmetry. Increase matrix column counter
+        // accordingly
+        s[1] = 0;
+      }
+    }
+    
+    i++;
+    
+  }
+  mdTi->print();
+  // md container
+  //md.clear();// this is to be used only for 1st order!!
+  md.remove(Symmetry(ord));
+  md.insert(mdTi);//(&mdTi);
+}
+
+
+void KordpDynare::writeModelInfo(Journal& jr) const
+{
+  // write info on variables
+  {
+    JournalRecordPair rp(journal);
+    rp << "Information on variables" << endrec;
+    JournalRecord rec1(journal);
+    rec1 << "Number of endogenous:            " << ny() << endrec;
+    JournalRecord rec2(journal);
+    rec2 << "Number of exogenous:             " << nexog() << endrec;
+    JournalRecord rec3(journal);
+    rec3 << "Number of static:                " << nstat() << endrec;
+    JournalRecord rec4(journal);
+    rec4 << "Number of predetermined:         " << npred()+nboth() << endrec;
+    JournalRecord rec5(journal);
+    rec5 << "Number of forward looking:       " << nforw()+nboth() << endrec;
+    JournalRecord rec6(journal);
+    rec6 << "Number of both:                  " << nboth() << endrec;
+  }
+  
+}
+/*********************************************************
+* LLxSteady()
+* returns ySteady extended with leads and lags suitable for 
+* passing to <model>_dynamic DLL
+*************************************************************/
+Vector * KordpDynare::LLxSteady( const Vector& yS){
+  if ((nJcols-nExog) == yS.length()) {
+    mexPrintf("k_ord_dynare.cpp: Warning in LLxSteady: ySteady already. right size");
+    return NULL;
+  }
+  // create temporary square 2D matrix size nEndo x nEndo (sparse)  
+  // for the lag, current and lead blocks of the jacobian
+  Vector * llxSteady = new Vector(nJcols-nExog); 
+  try{
+    for (int ll_row=0; ll_row< ll_Incidence->nrows(); ll_row++)
+    {
+      // populate (non-sparse) vector with ysteady values 
+      for (int i=0;i<nY;i++){
+        if(ll_Incidence->get(ll_row,i))
+          (*llxSteady)[((int)ll_Incidence->get(ll_row,i))-1] = yS[i];
+      }
+    }
+  } catch (const TLException& e) {
+    mexPrintf("Caugth TL exception in LLxSteady: ");
+    e.print();
+    return NULL;// 255;
+  }catch (...){
+    mexPrintf(" Error in LLxSteady - wrong index?");
+  }
+  
+#ifdef DEBUG		
+  for (int j=0;j<nJcols-nExog;j++)
+    mexPrintf("LLxSteady: [%d] =%f .\n", j, (*llxSteady)[j]);
+#endif	
+  return llxSteady;
+}
+
+
+/************************************
+* Reorder DynareJacobianIndices of variables in a vector according to
+* given int * varOrder together with lead & lag incidence matrix and 
+* any the extra columns for exogenous vars, and then, 
+* reorders its blocks given by the varOrder and the Dynare++ expectations:
+
+  * extra	nboth+ npred (t-1) lags
+  * varOrder
+		static:
+    pred
+    both    
+    forward    
+    * extra both + nforw (t+1) leads, and 
+    * extra exogen
+    
+      * so to match the jacobian organisation expected by the Appoximation class 
+      both + nforw (t+1) leads
+      static
+      pred
+      both
+      forward
+      nboth+ npred  (t-1) lags
+      exogen
+************************************/
+
+vector<int> * KordpDynare::ReorderDynareJacobianIndices( const vector<int> * varOrder){
+  // create temporary square 2D matrix size nEndo x nEndo (sparse)  
+  // for the lag, current and lead blocks of the jacobian
+  //	int * JacobianIndices = (int*) calloc(nJcols+1, sizeof(int)); 
+  vector<int> * JacobianIndices = new vector<int>(nJcols); 
+  vector <int> tmp(nY); 
+  int i,j, rjoff=nJcols-nExog-1; //, ll_off, j;
+#ifdef DEBUG		
+  mexPrintf("ReorderDynareJacobianIndice:ll_Incidence->nrows() =%d .\n", ll_Incidence->nrows());
+#endif
+  try{
+    for (int ll_row=0; ll_row< ll_Incidence->nrows(); ll_row++)
+    {
+      // reorder in orde-var order & populate temporary nEndo (sparse) vector with 
+      // the lag, current and lead blocks of the jacobian respectively
+      for (i=0;i<nY;i++){
+        //				tmp[varOrder[j]]=(int)ll_Incidence->get(ll_row,j); 
+        tmp[i]=((int)ll_Incidence->get(ll_row,(*varOrder)[i]-1)); 
+#ifdef DEBUG		
+        mexPrintf("get(ll_row,(*varOrder)[%d]-1)) = tmp[%d]=%d .\n", 
+          i, i, (int)ll_Incidence->get(ll_row,(*varOrder)[i]-1));
+#endif
+      }
+      // write the reordered blocks back to the jacobian
+      // in reverse order
+      for (j=nY-1;j>=0;j--){
+        if (tmp[j]){
+          (*JacobianIndices)[rjoff]=tmp[j] -1;
+          rjoff--;
+          if (rjoff<0){
+            //					mexPrintf(" Error in ReorderIndices - negative rjoff index?");
+            break;
+            //					return NULL;
+          }
+        }
+      }
+    }
+  } catch (const TLException& e) {
+    mexPrintf("Caugth TL exception in ReorderIndices: ");
+    e.print();
+    return NULL;// 255;
+  }catch (...){
+    mexPrintf(" Error in ReorderIndices - wrong index?");
+  }
+  //add the indices for the nExog exogenous jacobians
+  for (j=nJcols-nExog;j<nJcols;j++){
+    (*JacobianIndices)[j]=j;
+  }
+#ifdef DEBUG		
+  for (j=0;j<nJcols;j++)
+    mexPrintf("ReorderDynareJacobianIndice: [%d] =%d .\n", j, (*JacobianIndices)[j]);
+#endif	
+  return JacobianIndices;
+}
+
+
+/************************************
+* Reorder first set of columns of variables in a (jacobian) matrix 
+* according to order given in  varsOrder together with the extras 
+* assuming tdx ncols() - nExog is eaqual or less than length of varOrder and
+* of any of its elements too.
+************************************/
+
+void KordpDynare::ReorderCols(TwoDMatrix * tdx, const vector<int> * vOrder){
+  
+  if (tdx->ncols() > vOrder->size()){
+    mexPrintf(" Error in ReorderColumns - size of order var is too small");
+    return;
+  }
+  TwoDMatrix tmp(*tdx); // temporary 2D matrix
+  TwoDMatrix &tmpR=tmp;
+  tdx->zeros();// empty original matrix
+  // reorder the columns
+  try{
+    for (int i =0; i<tdx->ncols() ; i++)
+      tdx->copyColumn(tmpR,(*vOrder)[i],i);
+  } catch (const TLException& e) {
+    printf("Caugth TL exception in ReorderColumns: ");
+    e.print();
+    return;// 255;
+  }catch (...){
+    mexPrintf(" Error in ReorderColumns - wrong index?");
+  }
+}
+void KordpDynare::ReorderCols(TwoDMatrix * tdx, const int * vOrder){
+  
+  TwoDMatrix tmp(*tdx); // temporary 2D matrix
+  TwoDMatrix &tmpR=tmp;
+  tdx->zeros();// empty original matrix
+  // reorder the columns
+  try{
+    for (int i =0; i<tdx->ncols() ; i++)
+      tdx->copyColumn(tmpR,vOrder[i],i);
+  } catch (const TLException& e) {
+    printf("Caugth TL exception in ReorderColumns: ");
+    e.print();
+    return;// 255;
+  }catch (...){
+    mexPrintf(" Error in ReorderColumns - wrong index?");
+  }
+}
+
+/***********************************************************************
+* Recursive hierarchical block reordering of the higher order, input model 
+*	derivatives inc. Hessian   
+* This is now obsolete but kept in in case it is needed
+***********************************************************************/
+
+void KordpDynare::ReorderBlocks(TwoDMatrix * tdx, const vector<int> * vOrder){
+  // determine order of the matrix
+  
+  double dbOrder = log(tdx->ncols())/log(nJcols);
+  int ibOrder= (int) dbOrder;
+  if ((double )ibOrder != dbOrder || ibOrder>nOrder) {
+    mexPrintf(" Error in ReorderBlocks - wrong order %d", dbOrder);
+    return;
+  }
+  TwoDMatrix tmp(*tdx); // temporary 2D matrix
+  TwoDMatrix &tmpR=tmp;
+  tdx->zeros();// empty original matrix
+  
+  if(ibOrder>1){
+    int nBlocks=tmp.ncols()/ nJcols; //pow((float)nJcols,ibOrder-1);
+    int bSize=tmp.ncols()/nBlocks;
+    for (int j = 0; j<nBlocks;  ++j){
+      TwoDMatrix subtdx(tmpR, bSize*((*vOrder)[j]), bSize);
+      ReorderBlocks(&subtdx, vOrder);
+      tdx->place(subtdx, 0, bSize*j);
+    }
+  } else{
+    //ReorderColumns(TwoDMatrix * subtdx, const vector<int> * vOrder)
+    if (tdx->ncols() > vOrder->size()){
+      mexPrintf(" Error in ReorderColumns - size of order var is too small");
+      return;
+    }
+    // reorder the columns
+    try{
+      for (int i =0; i<tdx->ncols() ; i++)
+        tdx->copyColumn(tmpR,(*vOrder)[i],i);
+    } catch (const TLException& e) {
+      printf("Caugth TL exception in ReorderColumns: ");
+      e.print();
+      return;// 255;
+    }catch (...){
+      mexPrintf(" Error in ReorderColumns - wrong index?");
+    }
+  }
+}
+
+
+/****************************
+*  K-Order Perturbation instance of Jacobian:
+************************************/
+KordpJacobian::KordpJacobian(KordpDynare& dyn)
+: Jacobian(dyn.ny()), dyn(dyn)
+{
+  zeros();
+};
+
+void KordpJacobian::eval(const Vector& yy)
+{
+		dyn.calcDerivatives( yy, *this);
+    
+};
+
+void KordpVectorFunction::eval(const ConstVector& in, Vector& out)
+{
+  check_for_eval(in, out);
+  Vector xx(d.nexog());
+  xx.zeros();
+  d.evaluateSystem(out, in, xx);
+}
+
+/**************************************************************************************/
+/*       DynareNameList class                                                         */
+/**************************************************************************************/
+vector<int> DynareNameList::selectIndices(const vector<const char*>& ns) const
+{
+  vector<int> res;
+  for (unsigned int i = 0; i < ns.size(); i++) {
+    int j = 0;
+    while (j < getNum() && strcmp(getName(j), ns[i]) != 0)
+      j++;
+    if (j == getNum())
+      throw DynareException(__FILE__, __LINE__,
+      string("Couldn't find name for ") + ns[i] +
+      " in DynareNameList::selectIndices");
+    res.push_back(j);
+  }
+  return res;
+}
+
+DynareNameList::DynareNameList(const  KordpDynare& dynare)
+{
+  for (int i = 0; i < dynare.ny(); i++) {
+    names.push_back(dynare.dnl->getName(i));
+  }
+}
+DynareNameList::DynareNameList(const KordpDynare& dynare, const char ** namesp)
+{
+#ifdef DEBUG		
+  mexPrintf("k_ord_dynare DynareNameList.\n");
+  mexPrintf("k_ord_dynare DynareNameList dynare.ny=%d .\n", dynare.ny());
+#endif		
+  
+  for (int i = 0; i < dynare.ny(); i++) {
+#ifdef DEBUG		
+    mexPrintf("k_ord_dynare DynareNameList names[%d]= %s.\n", i, namesp[i] );
+#endif		
+    names.push_back(namesp[i]);
+  }
+}
+
+DynareExogNameList::DynareExogNameList(const KordpDynare& dynare)
+{
+  for (int i = 0; i < dynare.nexog(); i++) {
+    names.push_back(dynare.denl->getName(i));
+  }
+}
+
+DynareExogNameList::DynareExogNameList(const KordpDynare& dynare, const char ** namesp)
+{
+#ifdef DEBUG		
+  mexPrintf("k_ord_dynare DynareExogNameList dynare.nexog=%d .\n", dynare.nexog());
+#endif		
+  for (int i = 0; i < dynare.nexog(); i++) {
+#ifdef DEBUG		
+    mexPrintf("k_ord_dynare DynareExogNameList names[%d]= %s.\n", i, namesp[i] );
+#endif		
+    names.push_back(namesp[i]);
+  }
+}
+
+DynareStateNameList::DynareStateNameList(const KordpDynare& dynare, const DynareNameList& dnl,
+                                         const DynareExogNameList& denl)
+{
+  for (int i = 0; i < dynare.nys(); i++){
+#ifdef DEBUG		
+    mexPrintf("k_ord_dynare DynareStateNameList dnl names[%d]= %s.\n", i, dnl.getName(i+dynare.nstat()) );
+#endif		
+    names.push_back(dnl.getName(i+dynare.nstat()));
+  }
+  for (int i = 0; i < dynare.nexog(); i++){
+#ifdef DEBUG		
+    mexPrintf("k_ord_dynare DynareStateNameList denl names[%d]= %s.\n", i, denl.getName(i));
+#endif		
+    names.push_back(denl.getName(i));
+  }
+}
+
diff --git a/mex/sources/korderpert/src/k_ord_dynare.h b/mex/sources/korderpert/src/k_ord_dynare.h
index b1be48b3b..821a146e6 100644
--- a/mex/sources/korderpert/src/k_ord_dynare.h
+++ b/mex/sources/korderpert/src/k_ord_dynare.h
@@ -1,283 +1,250 @@
-/*
- * Copyright (C) 2005 Ondra Kamenik
- * Copyright (C) 2008-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 <http://www.gnu.org/licenses/>.
- */
-
-// based on: $Id: dynare3.h 1764 2008-03-31 14:30:55Z kamenik $
-// by 2005, Ondra Kamenik
-
-#ifndef K_ORD_DYNARE3_H
-#define K_ORD_DYNARE3_H
-#include <vector>
-#include "t_container.h"
-#include "sparse_tensor.h"
-#include "decision_rule.h"
-#include "dynamic_model.h"
-
-#include "exception.h"
-#include "dynare_exception.h"
-#include "fs_tensor.h"
-#include "SylvException.h"
-#include "tl_exception.h"
-#include "kord_exception.h"
-#include "nlsolve.h"
-#include "approximation.h"
-#include "k_order_perturbation.h"
-
-class KordpDynare;
-
-// derive from Approximation to get protected derivatives
-/******
- class FistOrderApproximation: public Approximation{
-    TwoDMatrix *gy;
-    TwoDMatrix *gu;
-    public:
-	FistOrderApproximation();
-	FistOrderApproximation(FistOrderApproximation& fo): ::Approximation(fo){
-		if (&(fo.GetGy())!=0){
-			gy= new TwoDMatrix(fo.GetGy());
-			gu= new TwoDMatrix(fo.GetGu());
-		}
-	};
-	virtual ~FistOrderApproximation(){
-		delete gy;
-		delete gu;
-	};
-	virtual void approxAtSteady(); // NOTE: change the Approximation parent to use virual too so that it can be overriden !!
-	//FGSContainer* GetRuleDers(){return rule_ders;};
-	//FGSContainer* GetRuleDersSS(){return rule_ders_ss;};
-	const TwoDMatrix& GetGy(){return (const TwoDMatrix&)(*gy);};
-	const TwoDMatrix& GetGu(){return (const TwoDMatrix&)(*gu);};
-    virtual void saveRuleDerivs(FirstOrder& fo); 
-};
-***********/
-
-
-
-/*////////////////////////////////////////////*/
-// instantiations of pure abstract class NameList in dynamic_model.h:
-/*////////////////////////////////////////////*/
-class DynareNameList : public NameList {
-    vector<const char*> names;
-public:
-    DynareNameList(const KordpDynare& dynare);
-    DynareNameList(const KordpDynare& dynare, const char ** names);
-    int getNum() const {return (int)names.size();}
-    const char* getName(int i) const {return names[i];}
-/** This for each string of the input vector calculates its index
-* in the names. And returns the resulting vector of indices. If
-* the name cannot be found, then an exception is raised. */
-vector<int> selectIndices(const vector<const char*>& ns) const;
-};
-
-class DynareExogNameList : public NameList {
-    vector<const char*> names;
-public:
-    DynareExogNameList(const  KordpDynare& dynare);
-    DynareExogNameList(const KordpDynare& dynare, const char ** names);
-    int getNum() const
-    {return (int)names.size();}
-    const char* getName(int i) const
-    {return names[i];}
-};
-		
-class DynareStateNameList : public NameList {
-    vector<const char*> names;
-public:
-    DynareStateNameList(const KordpDynare& dynare, const DynareNameList& dnl,
-    const DynareExogNameList& denl);
-    int getNum() const
-    {return (int)names.size();}
-    const char* getName(int i) const
-    {return names[i];}
-};
-/*********************************************/
-// The following only implements DynamicModel with help of ogdyn::DynareModel
-// instantiation of pure abstract DynamicModel decl. in dynamic_model.h
-//class DynamicModelDLL;
-class KordpJacobian;
-class KordpDynare : public DynamicModel {
-friend class DynareNameList;
-friend class DynareExogNameList;
-friend class DynareStateNameList;
-friend class KordpDynareJacobian;
-friend class DynamicModelDLL;
-	//////////
-	const int nStat;
-	const int nBoth;
-	const int nPred;
-	const int nForw;
-	const int nExog;
-	const int nPar;
-	const int nYs; // ={npred + nboth ; }
-	const int nYss; // nyss ={ nboth + nforw ; }
-	const int nY;  // = num_endo={ nstat + npred + nboth + nforw ; }
-	const int nJcols; // no of jacobian columns= nExog+nEndo+nsPred+nsForw
-	const int nSteps;
-	const int nOrder;
-	Journal& journal;
-	///	DynamicModel* model;
-	///const char* modName;
-	Vector* ySteady;
-	Vector* params;
-	TwoDMatrix* vCov;
-	TensorContainer<FSSparseTensor> md; // ModelDerivatives
-	DynareNameList* dnl;  
-	DynareExogNameList* denl;
-	DynareStateNameList* dsnl;
-	const double ss_tol;
-	const vector<int>* varOrder;
-	const TwoDMatrix * ll_Incidence;
-	double qz_criterium;
-	vector<int> * JacobianIndices;
-public:
-	KordpDynare(const char** endo, int num_endo,
-		const char** exo, int num_exo, int num_par, //const char** par,
-		Vector* ySteady, TwoDMatrix* vCov, Vector* params, int nstat,int nPred, 
-		int nforw, int nboth, const int nJcols, const int nSteps, const int ord, 	//const char* modName,
-		Journal& jr, DynamicModelDLL& dynamicDLL, double sstol, 
-		const vector<int>* varOrder, const TwoDMatrix * ll_Incidence, 
-		double qz_criterium );
-
-	/** Makes a deep copy of the object. */
-	KordpDynare(const KordpDynare& dyn);
-	virtual ~KordpDynare();
-	int nstat() const
-	{return nStat;}
-	int nboth() const
-	{return nBoth;}
-	int npred() const
-	{return nPred;}
-	int nforw() const
-	{return nForw;}
-	int nexog() const
-	{return nExog;}
-	int nys() const 
-	{return nYs;}
-	int nyss() const
-	{return nYss;}
-	int ny() const
-	{return nY;}
-	int steps() const
-	{return nSteps;}
-	int order() const
-	{return nOrder;}
-	const NameList& getAllEndoNames() const
-	{return *dnl;}
-	const NameList& getStateNames() const
-	{return *dsnl;}
-	const NameList& getExogNames() const
-	{return *denl;}
-	const TwoDMatrix& getVcov() const
-	{return *vCov;}
-	Vector& getParams()
-	{return *params;}
-	
-	const TensorContainer<FSSparseTensor>& getModelDerivatives() const
-	{return md;}
-	const Vector& getSteady() const
-	{return *ySteady;}
-	Vector& getSteady()
-	{return *ySteady;}
-	///	const ogdyn::DynareModel& getModel() const
-	///		{return *model;}
-	
-	// here is true public interface
-	void solveDeterministicSteady(Vector& steady);
-	void solveDeterministicSteady()
-	{solveDeterministicSteady(*ySteady);}
-	void evaluateSystem(Vector& out, const Vector& yy, const Vector& xx);
-	void evaluateSystem(Vector& out, const Vector& yym, const Vector& yy,
-		const Vector& yyp, const Vector& xx);
-	void calcDerivatives(const Vector& yy, const Vector& xx);
-	//void calcDerivatives(const Vector& yy, TwoDMatrix& jj);
-    void calcDerivatives(const Vector& yy, ogu::Jacobian& jacob);
-	void calcDerivativesAtSteady();
-	DynamicModelDLL& dynamicDLL;
-	///	void writeMat4(FILE* fd, const char* prefix) const;
-	///	void writeDump(const std::string& basename) const;
-	DynamicModel* clone() const
-	{return new KordpDynare(*this);}
-	void ReorderCols(TwoDMatrix * tdx, const int * varOrder);
-	void ReorderCols(TwoDMatrix * tdx, const vector<int> * varOrder);
-	Vector * LLxSteady( const Vector& yS); // returns ySteady extended with leads and lags 
-
-private:
-	void writeModelInfo(Journal& jr) const;
-	int * ReorderDynareJacobianIndices( const int * varOrder);
-	vector<int> * ReorderDynareJacobianIndices( const vector<int> * varOrder);
-	void ReorderBlocks(TwoDMatrix * tdx, const int * varOrder);
-	void ReorderBlocks(TwoDMatrix * tdx, const vector<int> * vOrder);
-  void populateDerivativesContainer(TwoDMatrix*g, int ord, const vector<int>* vOrder);
-};
-
-/****************************
-* ModelDerivativeContainer manages derivatives container
-************************************/
-
-class ModelDerivativeContainer //: public ogp::FormulaDerEvalLoader 
-{
-protected:
-//	const ogp::FineAtoms& atoms;
-	TensorContainer<FSSparseTensor>& md;
-public:
-	ModelDerivativeContainer(const KordpDynare& model, TensorContainer<FSSparseTensor>& mod_ders,
-						int order);
-	void load(int i, int iord, const int* vars, double res);
-};
-
-  /****************************
-  *  K-Order Perturbation instance of Jacobian:
-  ************************************/
-  
-  class KordpJacobian : public ogu::Jacobian ///, public ogp::FormulaDerEvalLoader 
-  {
-  protected:
-	  KordpDynare& dyn;
-  public:
-	  KordpJacobian(	KordpDynare& dyn);
-	  virtual ~KordpJacobian() {}
-	  // Load <mod>_dynamic.DLL 
-//	  void load(const char** modName);
-	  void eval(const Vector& in);
-  };
-  
-  
-  /****************************
-  *  K-Order Perturbation instance of VectorFunction:
-  ************************************/
-  
-  class 	KordpVectorFunction : public ogu::VectorFunction {
-  protected:
-	  KordpDynare& d;
-  public:
-	  KordpVectorFunction(	KordpDynare& dyn)
-		  : d(dyn) {}
-	  virtual ~KordpVectorFunction() {}
-	  int inDim() const
-	  {return d.ny();}
-	  int outDim() const
-	  {return d.ny();}
-	  void eval(const ConstVector& in, Vector& out);
-  };
-  
-#endif
-  
-  // Local Variables:
-  // mode:C++
-  // End:
+/*
+* Copyright (C) 2005 Ondra Kamenik
+* Copyright (C) 2008-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 <http://www.gnu.org/licenses/>.
+*/
+
+// based on: $Id: dynare3.h 1764 2008-03-31 14:30:55Z kamenik $
+// by 2005, Ondra Kamenik
+
+#ifndef K_ORD_DYNARE3_H
+#define K_ORD_DYNARE3_H
+#include <vector>
+#include "t_container.h"
+#include "sparse_tensor.h"
+#include "decision_rule.h"
+#include "dynamic_model.h"
+
+#include "exception.h"
+#include "dynare_exception.h"
+#include "fs_tensor.h"
+#include "SylvException.h"
+#include "tl_exception.h"
+#include "kord_exception.h"
+#include "nlsolve.h"
+#include "approximation.h"
+#include "k_order_perturbation.h"
+
+class KordpDynare;
+
+/*////////////////////////////////////////////*/
+// instantiations of pure abstract class NameList in dynamic_model.h:
+/*////////////////////////////////////////////*/
+class DynareNameList : public NameList {
+  vector<const char*> names;
+public:
+  DynareNameList(const KordpDynare& dynare);
+  DynareNameList(const KordpDynare& dynare, const char ** names);
+  int getNum() const {return (int)names.size();}
+  const char* getName(int i) const {return names[i];}
+  /** This for each string of the input vector calculates its index
+  * in the names. And returns the resulting vector of indices. If
+  * the name cannot be found, then an exception is raised. */
+  vector<int> selectIndices(const vector<const char*>& ns) const;
+};
+
+class DynareExogNameList : public NameList {
+  vector<const char*> names;
+public:
+  DynareExogNameList(const  KordpDynare& dynare);
+  DynareExogNameList(const KordpDynare& dynare, const char ** names);
+  int getNum() const
+  {return (int)names.size();}
+  const char* getName(int i) const
+  {return names[i];}
+};
+
+class DynareStateNameList : public NameList {
+  vector<const char*> names;
+public:
+  DynareStateNameList(const KordpDynare& dynare, const DynareNameList& dnl,
+    const DynareExogNameList& denl);
+  int getNum() const
+  {return (int)names.size();}
+  const char* getName(int i) const
+  {return names[i];}
+};
+/*********************************************/
+// The following only implements DynamicModel with help of ogdyn::DynareModel
+// instantiation of pure abstract DynamicModel decl. in dynamic_model.h
+//class DynamicModelDLL;
+class KordpJacobian;
+class KordpDynare : public DynamicModel {
+  friend class DynareNameList;
+  friend class DynareExogNameList;
+  friend class DynareStateNameList;
+  friend class KordpDynareJacobian;
+  friend class DynamicModelDLL;
+  //////////
+  const int nStat;
+  const int nBoth;
+  const int nPred;
+  const int nForw;
+  const int nExog;
+  const int nPar;
+  const int nYs; // ={npred + nboth ; }
+  const int nYss; // nyss ={ nboth + nforw ; }
+  const int nY;  // = num_endo={ nstat + npred + nboth + nforw ; }
+  const int nJcols; // no of jacobian columns= nExog+nEndo+nsPred+nsForw
+  const int nSteps;
+  const int nOrder;
+  Journal& journal;
+  ///	DynamicModel* model;
+  ///const char* modName;
+  Vector* ySteady;
+  Vector* params;
+  TwoDMatrix* vCov;
+  TensorContainer<FSSparseTensor> md; // ModelDerivatives
+  DynareNameList* dnl;  
+  DynareExogNameList* denl;
+  DynareStateNameList* dsnl;
+  const double ss_tol;
+  const vector<int>* varOrder;
+  const TwoDMatrix * ll_Incidence;
+  double qz_criterium;
+  vector<int> * JacobianIndices;
+public:
+  KordpDynare(const char** endo, int num_endo,
+    const char** exo, int num_exo, int num_par, //const char** par,
+    Vector* ySteady, TwoDMatrix* vCov, Vector* params, int nstat,int nPred, 
+    int nforw, int nboth, const int nJcols, const int nSteps, const int ord, 	//const char* modName,
+    Journal& jr, DynamicModelDLL& dynamicDLL, double sstol, 
+    const vector<int>* varOrder, const TwoDMatrix * ll_Incidence, 
+    double qz_criterium );
+  
+  /** Makes a deep copy of the object. */
+  KordpDynare(const KordpDynare& dyn);
+  virtual ~KordpDynare();
+  int nstat() const
+  {return nStat;}
+  int nboth() const
+  {return nBoth;}
+  int npred() const
+  {return nPred;}
+  int nforw() const
+  {return nForw;}
+  int nexog() const
+  {return nExog;}
+  int nys() const 
+  {return nYs;}
+  int nyss() const
+  {return nYss;}
+  int ny() const
+  {return nY;}
+  int steps() const
+  {return nSteps;}
+  int order() const
+  {return nOrder;}
+  const NameList& getAllEndoNames() const
+  {return *dnl;}
+  const NameList& getStateNames() const
+  {return *dsnl;}
+  const NameList& getExogNames() const
+  {return *denl;}
+  const TwoDMatrix& getVcov() const
+  {return *vCov;}
+  Vector& getParams()
+  {return *params;}
+  
+  const TensorContainer<FSSparseTensor>& getModelDerivatives() const
+  {return md;}
+  const Vector& getSteady() const
+  {return *ySteady;}
+  Vector& getSteady()
+  {return *ySteady;}
+  
+  // here is true public interface
+  void solveDeterministicSteady(Vector& steady);
+  void solveDeterministicSteady()
+  {solveDeterministicSteady(*ySteady);}
+  void evaluateSystem(Vector& out, const Vector& yy, const Vector& xx);
+  void evaluateSystem(Vector& out, const Vector& yym, const Vector& yy,
+    const Vector& yyp, const Vector& xx);
+  void calcDerivatives(const Vector& yy, const Vector& xx);
+  //void calcDerivatives(const Vector& yy, TwoDMatrix& jj);
+  void calcDerivatives(const Vector& yy, ogu::Jacobian& jacob);
+  void calcDerivativesAtSteady();
+  DynamicModelDLL& dynamicDLL;
+  ///	void writeMat4(FILE* fd, const char* prefix) const;
+  ///	void writeDump(const std::string& basename) const;
+  DynamicModel* clone() const
+  {return new KordpDynare(*this);}
+  void ReorderCols(TwoDMatrix * tdx, const int * varOrder);
+  void ReorderCols(TwoDMatrix * tdx, const vector<int> * varOrder);
+  Vector * LLxSteady( const Vector& yS); // returns ySteady extended with leads and lags 
+  
+private:
+  void writeModelInfo(Journal& jr) const;
+  int * ReorderDynareJacobianIndices( const int * varOrder);
+  vector<int> * ReorderDynareJacobianIndices( const vector<int> * varOrder);
+  void ReorderBlocks(TwoDMatrix * tdx, const int * varOrder);
+  void ReorderBlocks(TwoDMatrix * tdx, const vector<int> * vOrder);
+  void populateDerivativesContainer(TwoDMatrix*g, int ord, const vector<int>* vOrder);
+};
+
+/****************************
+* ModelDerivativeContainer manages derivatives container
+************************************/
+
+class ModelDerivativeContainer //: public ogp::FormulaDerEvalLoader 
+{
+protected:
+  //	const ogp::FineAtoms& atoms;
+  TensorContainer<FSSparseTensor>& md;
+public:
+  ModelDerivativeContainer(const KordpDynare& model, TensorContainer<FSSparseTensor>& mod_ders,
+    int order);
+  void load(int i, int iord, const int* vars, double res);
+};
+
+/****************************
+*  K-Order Perturbation instance of Jacobian:
+************************************/
+
+class KordpJacobian : public ogu::Jacobian ///, public ogp::FormulaDerEvalLoader 
+{
+protected:
+  KordpDynare& dyn;
+public:
+  KordpJacobian(	KordpDynare& dyn);
+  virtual ~KordpJacobian() {}
+  // Load <mod>_dynamic.DLL 
+  //	  void load(const char** modName);
+  void eval(const Vector& in);
+};
+
+
+/****************************
+*  K-Order Perturbation instance of VectorFunction:
+************************************/
+
+class 	KordpVectorFunction : public ogu::VectorFunction {
+protected:
+  KordpDynare& d;
+public:
+  KordpVectorFunction(	KordpDynare& dyn)
+		  : d(dyn) {}
+  virtual ~KordpVectorFunction() {}
+  int inDim() const
+  {return d.ny();}
+  int outDim() const
+  {return d.ny();}
+  void eval(const ConstVector& in, Vector& out);
+};
+
+#endif
+
diff --git a/mex/sources/korderpert/src/k_order_perturbation.cpp b/mex/sources/korderpert/src/k_order_perturbation.cpp
index a18e08e13..a6b6cf62a 100644
--- a/mex/sources/korderpert/src/k_order_perturbation.cpp
+++ b/mex/sources/korderpert/src/k_order_perturbation.cpp
@@ -1,625 +1,625 @@
-/*
- * Copyright (C) 2008-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 <http://www.gnu.org/licenses/>.
- */
-
-/****************************************************** 
-// k_order_perturbation.cpp : Defines the entry point for the k-order perturbation application DLL.
-//
-// called from Dynare dr1_k_order.m, (itself called form resol.m instead of regular dr1.m)
-//            if options_.order < 2 % 1st order only
-//                [ysteady, ghx_u]=k_order_perturbation(dr,task,M_,options_, oo_ , ['.' mexext]);
-//            else % 2nd order
-//                [ysteady, ghx_u, g_2]=k_order_perturbation(dr,task,M_,options_, oo_ , ['.' mexext]);
-// inputs: 
-//			dr,		- Dynare structure
-//			task,  - check or not, not used
-//			M_		- Dynare structure
-//			options_ - Dynare structure
-//			oo_		- Dynare structure 
-//			['.' mexext] Matlab dll extension
-// returns:
-//			 ysteady steady state
-//			ghx_u - first order rules packed in one matrix
-//			g_2 - 2nd order rules packed in one matrix
-**********************************************************/
-
-
-#include "k_ord_dynare.h"
-#include "math.h"
-#include <cstring>
-
-#include <cctype>
-
-#ifdef  _MSC_VER  //&&WINDOWS 
-
-
-BOOL APIENTRY DllMain( HANDLE hModule, 
-					  DWORD  ul_reason_for_call, 
-					  LPVOID lpReserved
-					  )
-{
-    switch (ul_reason_for_call)
-	{
-	case DLL_PROCESS_ATTACH:
-	case DLL_THREAD_ATTACH:
-	case DLL_THREAD_DETACH:
-	case DLL_PROCESS_DETACH:
-		break;
-    }
-    return TRUE;
-}
-
-// Some MS Windows preambles
-// This is an example of an exported variable
-K_ORDER_PERTURBATION_API int nK_order_perturbation=0;
-
-// This is an example of an exported function.
-K_ORDER_PERTURBATION_API int fnK_order_perturbation(void)
-{
-	return 42;
-}
-// This is the constructor of a class that has been exported.
-// see k_order_perturbation.h for the class definition
-CK_order_perturbation::CK_order_perturbation()
-{ 
-	return; 
-}
-
-#endif // _MSC_VER && WINDOWS
-
-extern "C" {
-	
-	// mexFunction: Matlab Inerface point and the main application driver	
-
-	void mexFunction(int nlhs, mxArray* plhs[],
-		int nrhs, const mxArray* prhs[])
-	{
-		if (nrhs < 5)
-			mexErrMsgTxt("Must have at least 5 input parameters.\n");
-		if (nlhs == 0)
-			mexErrMsgTxt("Must have at least 1 output parameter.\n");
-		
-		const mxArray* dr = prhs[0];
-		const int check_flag = (int)mxGetScalar(prhs[1]);
-		const mxArray* M_ = prhs[2];
-		const mxArray* options_=  (prhs[3]);
-		const mxArray* oo_ =  (prhs[4]);
-		
-		mxArray* mFname =mxGetField(M_, 0, "fname");
-		if(!mxIsChar(mFname)){
-			mexErrMsgTxt("Input must be of type char.");
-		}
-		const char*  fName = mxArrayToString(mFname);
-		const char * dfExt=NULL ;//Dyanamic file extension, e.g.".dll" or .mexw32; 
-		if (prhs[5] != NULL){
-			const mxArray* mexExt = prhs[5];
-			dfExt= mxArrayToString(mexExt);
-		} 
-
-#ifdef DEBUG		
-		mexPrintf("k_order_perturbation: check_flag = %d ,  fName = %s , mexExt=%s.\n", check_flag,fName, dfExt);
-#endif		
-        int kOrder;
-		mxArray* mxFldp = mxGetField(options_, 0,"order" );
-		if (mxIsNumeric(mxFldp))
-			kOrder = (int)mxGetScalar(mxFldp);
-		else
-			kOrder = 1;
-
-		double qz_criterium = 1+1e-6;
-		mxFldp = mxGetField(options_, 0,"qz_criterium" );
-		if (mxIsNumeric(mxFldp))
-			qz_criterium = (double)mxGetScalar(mxFldp);
-		
-		mxFldp 	= mxGetField(M_, 0,"params" );
-		double * dparams = (double *) mxGetData(mxFldp);
-		int npar = (int)mxGetM(mxFldp);
-		Vector * modParams =  new Vector(dparams, npar);
-#ifdef DEBUG		
-    mexPrintf("k_ord_perturbation: qz_criterium=%g, nParams=%d .\n",qz_criterium,npar);  
-    for (int i = 0; i < npar; i++) {
-        mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]);  }
-//	for (int i = 0; i < npar; i++) {
-//        mexPrintf("k_ord_perturbation: params_vec[%d]= %g.\n", i, params_vec[i] );   }
-#endif		        
-
-		mxFldp 	= mxGetField(M_, 0,"Sigma_e" );
-		dparams = (double *) mxGetData(mxFldp);
-		npar = (int)mxGetN(mxFldp);
-		TwoDMatrix * vCov =  new TwoDMatrix(npar, npar, dparams);
-
-
-//		mxFldp 	= mxGetField(oo_, 0,"steady_state" ); // use in order of declaration
-		mxFldp 	= mxGetField(dr, 0,"ys" );  // and not in order of dr.order_var
-//		mxFldp 	= mxGetField(oo_, 0,"dyn_ys" );  // and NOT extended ys
-		dparams = (double *) mxGetData(mxFldp);
-		const int nSteady = (int)mxGetM(mxFldp);
-		Vector * ySteady =  new Vector(dparams, nSteady);
-
-
-		mxFldp = mxGetField(dr, 0,"nstatic" );
-		const int nStat = (int)mxGetScalar(mxFldp);
-		mxFldp = mxGetField(dr, 0,"npred" );
-		int nPred = (int)mxGetScalar(mxFldp);
-		mxFldp = mxGetField(dr, 0,"nspred" );
-		const int nsPred = (int)mxGetScalar(mxFldp);
-		mxFldp = mxGetField(dr, 0,"nboth" );
-		const int nBoth = (int)mxGetScalar(mxFldp);
-		mxFldp = mxGetField(dr, 0,"nfwrd" );
-		const int nForw = (int)mxGetScalar(mxFldp);
-		mxFldp = mxGetField(dr, 0,"nsfwrd" );
-		const int nsForw = (int)mxGetScalar(mxFldp);
-
-		mxFldp = mxGetField(M_, 0,"exo_nbr" );
-		const int nExog = (int)mxGetScalar(mxFldp);
-		mxFldp = mxGetField(M_, 0,"endo_nbr" );
-		const int nEndo = (int)mxGetScalar(mxFldp);
-		mxFldp = mxGetField(M_, 0,"param_nbr" );
-		const int nPar = (int)mxGetScalar(mxFldp);
-        // it_ should be set to M_.maximum_lag
-		mxFldp = mxGetField(M_, 0,"maximum_lag" );
-		const int nMax_lag = (int)mxGetScalar(mxFldp);
-
-		nPred -= nBoth; // correct nPred for nBoth.
-
-		mxFldp 	= mxGetField(dr, 0,"order_var" );
-		dparams = (double *) mxGetData(mxFldp);
-		npar = (int)mxGetM(mxFldp);
-		if (npar != nEndo) {  //(nPar != npar)
-			mexErrMsgTxt("Incorrect number of input var_order vars.\n");
-			//return;
-		} 
-		vector<int> * var_order_vp = (new vector<int>(nEndo));//nEndo));
-		for (int v=0;v<nEndo;v++){
-			(*var_order_vp)[v] =(int)(*(dparams++));//[v];
-#ifdef DEBUG		
-			mexPrintf("var_order_vp)[%d] = %d .\n", v, (*var_order_vp)[v] );
-#endif
-		}
-
-		// the lag, current and lead blocks of the jacobian respectively
-		mxFldp 	= mxGetField(M_, 0,"lead_lag_incidence" );
-		dparams = (double *) mxGetData(mxFldp);
-		npar = (int)mxGetN(mxFldp);
-		int nrows = (int)mxGetM(mxFldp);
-#ifdef DEBUG		
-		mexPrintf("ll_Incidence nrows=%d, ncols = %d .\n", nrows, npar);
-#endif
-		TwoDMatrix * llincidence =  new TwoDMatrix ( nrows, npar,dparams);
-		if (npar != nEndo){
-			mexPrintf("Incorrect length of lead lag incidences: ncol=%d != nEndo =%d .\n", npar, nEndo);
-			return;
-		}
-#ifdef DEBUG		
-		for (int j=0;j<nrows;j++){
-			for (int i=0;i<nEndo;i++){
-				mexPrintf("llincidence->get(%d,%d) =%d .\n", 
-					j, i, (int)llincidence->get(j,i));}}
-#endif
-
-        const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
-        mexPrintf("k_order_perturbation: jcols= %d .\n", jcols);
-
-        mxFldp= mxGetField(M_, 0,"var_order_endo_names" );
-        mexPrintf("k_order_perturbation: Get nendo .\n");
-        const int nendo = (int)mxGetM(mxFldp);
-        const int widthEndo = (int)mxGetN(mxFldp);
-        const char**  endoNamesMX= DynareMxArrayToString( mxFldp,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 = (int)mxGetM(mxFldp);
-		const int widthExog = (int)mxGetN(mxFldp);
-        const char**  exoNamesMX= DynareMxArrayToString( mxFldp,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)
-			mexErrMsgTxt("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[%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 = 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: Calling dynamicDLL constructor.\n");
-#endif				
-			//			DynamicFn * pDynamicFn = loadModelDynamicDLL (fname);
-			DynamicModelDLL dynamicDLL(fName,nEndo, jcols, nMax_lag, nExog, dfExt);
-
-			// intiate tensor library
-#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 dynare constructor .\n");
-#endif			
-			// make KordpDynare object
-			KordpDynare dynare(endoNamesMX,  nEndo, exoNamesMX,  nExog, nPar, // paramNames,
-   			   ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,
-			   jcols, nSteps, kOrder, journal, dynamicDLL, sstol, var_order_vp, 
-			   llincidence,qz_criterium );
-
-			// construct main K-order approximation class
-#ifdef DEBUG		
-	mexPrintf("k_order_perturbation: Call Approximation constructor with qz_criterium=%f \n", qz_criterium);
-#endif
-			Approximation app(dynare, journal,  nSteps, false, qz_criterium);
-            // run stochastic steady 
-#ifdef DEBUG		
-	mexPrintf("k_order_perturbation: Calling walkStochSteady.\n");
-#endif			
-            app.walkStochSteady();		
-
-
-			ConstTwoDMatrix ss(app.getSS());
-			// 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);
-			}
-
-			std::string ss_matrix_name(fName);//params.prefix);
-			ss_matrix_name += "_steady_states";
-//			ConstTwoDMatrix(app.getSS()).writeMat4(matfd, ss_matrix_name.c_str());
-			ss.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);
-
-			/* Write derivative outputs into memory map */
-			map<string,ConstTwoDMatrix> mm;
-			app.getFoldDecisionRule().writeMMap(&mm);
-
-#ifdef DEBUG		
-			app.getFoldDecisionRule().print();
-			mexPrintf("k_order_perturbation: Map print: \n");
-			for (map<string,ConstTwoDMatrix>::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());
-            
-
-		// developement of the output.
-#ifdef DEBUG		
-		mexPrintf("k_order_perturbation: Filling 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. */
-				TwoDMatrix tmp_ss(nEndo,1, mxGetPr(plhs[0]));
-				tmp_ss = (const TwoDMatrix&)ss;
-#ifdef DEBUG		
-//				tmp_ss.print();  !! This print Crashes???
-#endif					
-			}
-			if (nlhs >= 2)
-			{
-				/* Set the output pointer to the combined output matrix gyu = [gy gu]. */
-				int ii=1;
-				for (map<string,ConstTwoDMatrix>::const_iterator cit=mm.begin();
-					((cit !=mm.end())&&(ii<nlhs)); ++cit) 
-				{
-					//if ((*cit).first!="g_0") 
-					{
-						plhs[ii] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL);
-						TwoDMatrix dgyu((*cit).second.numRows(), (*cit).second.numCols(), mxGetPr(plhs[ii]));
-						dgyu = (const TwoDMatrix&)(*cit).second;
-#ifdef DEBUG		
-			mexPrintf("k_order_perturbation: cit %d print: \n", ii);
-                        (*cit).second.print();
-			mexPrintf("k_order_perturbation: dguy %d print: \n", ii);
-  //                      dgyu.print(); !! This print Crashes???
-#endif					
-						++ii;
-					}
-				}
-			}
-		} catch (const KordException& e) {
-			printf("Caugth Kord exception: ");
-			e.print();
-			mexPrintf("Caugth Kord exception: %s", e.get_message());
-			std::string errfile(fName);//(params.basename);
-			errfile += "_error.log";
-			FILE* errfd = NULL;
-			if (NULL == (errfd=fopen(errfile.c_str(), "wb"))) {
-				fprintf(stderr, "Couldn't open %s for writing.\n", errfile.c_str());
-				return;// e.code();
-			}
-			fprintf(errfd ,"Caugth Kord exception: %s", e.get_message());
-			fclose(errfd);
-			return;// e.code();
-		} catch (const TLException& e) {
-			printf("Caugth TL exception: ");
-			e.print();
-			return;// 255;
-		} catch (SylvException& e) {
-			printf("Caught Sylv exception: ");
-			e.printMessage();
-			return;// 255;
-		} catch (const DynareException& e) {
-			printf("Caught KordpDynare exception: %s\n", e.message());
-			mexPrintf("Caugth Dynare exception: %s", e.message());
-			std::string errfile(fName);//(params.basename);
-			errfile += "_error.log";
-			FILE* errfd = NULL;
-			if (NULL == (errfd=fopen(errfile.c_str(), "wb"))) {
-				fprintf(stderr, "Couldn't open %s for writing.\n", errfile.c_str());
-				return;// e.code();
-			}
-			fprintf(errfd ,"Caugth KordDynare  exception: %s", e.message());
-			fclose(errfd);
-			return;// 255;
-		} catch (const ogu::Exception& e) {
-			printf("Caught ogu::Exception: ");
-			e.print();
-			mexPrintf("Caugth general exception: %s", e.message());
-			std::string errfile(fName);//(params.basename);
-			errfile += "_error.log";
-			FILE* errfd = NULL;
-			if (NULL == (errfd=fopen(errfile.c_str(), "wb"))) {
-				fprintf(stderr, "Couldn't open %s for writing.\n", errfile.c_str());
-				return;// e.code();
-			}
-			e.print(errfd);
-			fclose(errfd);
-			return;// 255;
-        }
-    };
-}; // end of extern C
-
-//////////////////////////////////////////////////////
-// Convert Matlab Dynare endo and exo names array to C type array of string pointers
-// Poblem is that Matlab mx function returns a long string concatenated by columns rather than rows
-// hence a rather low level approach is needed
-///////////////////////////////////////////////////////
-const char ** DynareMxArrayToString(const mxArray * mxFldp,const int len,const int width )
-{
-            char * cNamesCharStr= mxArrayToString(mxFldp);
-			const char ** ret = DynareMxArrayToString(cNamesCharStr,len, width );
-			return ret;
-}
-
-const char ** DynareMxArrayToString(const char * cNamesCharStr,const int len,const int width )
-{
-             char cNamesMX[len][width+1] ;//
-#ifdef DEBUG
-mexPrintf("loop DynareMxArrayToString cNamesCharStr = %s \n", cNamesCharStr);
-#endif
-            for (int i=0;i<width;i++){
-                for (int j=0;j<len;j++){
-     			// Allow alphanumeric and underscores "_" only:
-                    if (isalnum(cNamesCharStr[j+i*len])||('_'==cNamesCharStr[j+i*len])){
-                        cNamesMX[j][i]=cNamesCharStr[j+i*len];
-                    }
-                    else cNamesMX[j][i]='\0';
-                }
-            }
-            const char ** ret= (const char **)mxCalloc (len, sizeof(char*));
-            for (int j=0;j<len;j++){
-				cNamesMX[j][width]='\0';
-#ifdef DEBUG
-//				mexPrintf("String [%d]= %s \n", j, cNamesMX[j]);
-#endif
-                char * token= (char*) mxCalloc ( strlen(cNamesMX[j])+1,sizeof(char));
-                strcpy (token, cNamesMX[j]);
-				ret[j]=token;
-#ifdef DEBUG
-				mexPrintf("ret [%d]= %s \n", j, ret[j]);
-#endif
-			}
-			return ret;
-}
-
-
-/***********************************
-* Members of DynamicModelDLL for handling loading and calling 
-* <model>_dynamic () function
-**************************************/
-DynamicModelDLL::DynamicModelDLL(const char * modName, const int y_length, const int j_cols, 
-								 const int n_max_lag, const int n_exog, const char * sExt)
-	: length(y_length),jcols( j_cols), nMax_lag(n_max_lag), nExog(n_exog)
-{
-    char fName[MAX_MODEL_NAME];
-    strcpy(fName,modName);
-	using namespace std;
-	strcat(fName,"_dynamic");
-#ifdef DEBUG
-		mexPrintf("MexPrintf: Call Load run DLL %s .\n", fName);
-#endif	
-
-	try {
-#ifdef WINDOWS
-		if (sExt==NULL) sExt=(".dll");
-		HINSTANCE dynamicHinstance;
-//		dynamicHinstance=::LoadLibraryEx(strcat(fNname,"_.dll"),NULL,DONT_RESOLVE_DLL_REFERENCES);//sExt); //"_.dll");
-		dynamicHinstance=::LoadLibrary(strcat(fName,sExt));//.dll); //"_.dll");
-		if (dynamicHinstance==NULL)
-			throw 1; //alt: return;
-		//		(DynamicFn*)	typedef void * (__stdcall *DynamicFn)();
-#ifdef DEBUG
-		mexPrintf("MexPrintf: Call GetProcAddress  %s .\n", fName);
-#endif
-		Dynamic = (DynamicFn *) ::GetProcAddress(dynamicHinstance,"Dynamic");
-        
-# else // __linux__
-		if (sExt==NULL) sExt=(".so");
-		dynamicHinstance = dlopen(strcat(fName,sExt), RTLD_NOW);
-		if((dynamicHinstance == NULL) || dlerror()){
-			cerr << dlerror() << endl;
-			mexPrintf("MexPrintf:Error loading DLL: %s", dlerror);
-			throw 1;
-		}
-		Dynamic = (DynamicFn)dlsym(dynamicHinstance, "Dynamic");
-		if((Dynamic  == NULL) || dlerror()){
-			cerr << dlerror() << endl;
-			mexPrintf("MexPrintf:Error finding DLL function: %s", dlerror);
-			throw 2;
-		}
-# endif
-		
-	} catch (int i) {
-		mexPrintf("MexPrintf: error in Load and run DLL %s , %d.\n", fName, i);
-		mexErrMsgTxt("Err: An error in Load and run DLL  .\n");
-		return;
-		
-	} catch (...) {
-		mexPrintf("MexPrintf: Unknown error in Call MATLAB function %s.\n", fName);
-		mexErrMsgTxt("Err: Unknown error in Load and run DLL  .\n");
-		return;
-	}
-}
-
-// close DLL: If the referenced object was successfully closed, 
-// close() returns 0, non 0 otherwise
-int DynamicModelDLL::close(){
-#ifdef WINDOWS
-    // MS FreeLibrary returns non 0 if OK, 0 if fails.
-	bool rb=FreeLibrary(dynamicHinstance);
-	if (rb)
-		return 0;
-	else 
-		return 1;
-# else // linux
-	//If OK, dlclose() returns 0, non 0 otherwise
-	return dlclose(dynamicHinstance);
-# endif
-};
-
-
-void DynamicModelDLL::eval(const Vector&y, const TwoDMatrix&x, const  Vector* modParams, 
-		int it_, Vector&residual, TwoDMatrix*g1, TwoDMatrix*g2){
-
-        double  *dresidual, *dg1=NULL, *dg2=NULL; 
-        //int length=y.length(); // not!
-		if ((jcols-nExog)!=y.length()){
-			// throw DLL Error
-			mexPrintf(" DLL Error: (jcols-nExog)!=ys.length() \n");
-			return;
-		}
-        if (residual.length()<length){ // dummy or insufficient
-            Vector*tempv= new Vector(length );
-			residual=*tempv;
-            delete tempv;
-			residual.zeros();
-		}
-        if (g1!=NULL){
-            if (g1->nrows()!=length){ // dummy
-                delete g1;
-                g1=	new TwoDMatrix( length, jcols); // and get a new one
-				g1->zeros();
-            }
-            dg1= const_cast<double*>(g1->base());
-        }
-        if (g2!=NULL){
-            if (g2->nrows()!=length){ // dummy 
-                delete g2;
-                g2=	new TwoDMatrix( length, jcols*jcols);// and get a new one
-				g2->zeros();
-            }
-            dg2= const_cast<double*>(g2->base());
-        }
-        dresidual=const_cast<double*>(residual.base());
-        double *dy=const_cast<double*>(y.base());
-        double *dx=const_cast<double*>(x.base());
-        double *dbParams=const_cast<double*>(modParams->base());
-#ifdef DEBUG		
-    	mexPrintf(" try eval Dynamic with ne g1: cols=%d , rows=%d\n"
-            , g1->ncols(),g1->nrows());
-        for (int i = 0; i < modParams->length(); i++) {
-            mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]);  }
-        for (int i = 0; i < jcols-nExog; i++) {
-            mexPrintf("k_ord_perturbation: Ys[%d]= %g.\n", i, dy[i]);} 
-	mexPrintf("k_order_perturbation: call <model> Dynamic dParams= %g ,  , dy = %g dx = %f .\n"
-            ,dbParams[0],dy[0],dx[0]);
-
-#endif        
-        try{
-            Dynamic(dy, dx, nExog, dbParams, it_, dresidual, dg1, dg2);
-        }catch (...){
-            mexPrintf("MexPrintf: error in run Dynamic DLL \n");
-        }
-};
-
-void DynamicModelDLL::eval(const Vector&y, const TwoDMatrix&x, const Vector * modParams, 
-		Vector&residual, TwoDMatrix*g1, TwoDMatrix*g2){
-
-		eval(y, x, modParams, nMax_lag, residual, g1, g2);	
-};
-
-void DynamicModelDLL::eval(const Vector&y, const Vector&x, const Vector * modParams, 
-		Vector&residual, TwoDMatrix*g1, TwoDMatrix*g2){
-
-		/** ignore given exogens and create new 2D x matrix since 
-		* when calling <model>_dynamic(z,x,params,it_) x must be equal to
-		* zeros(M_.maximum_lag+1,M_.exo_nbr)
-		**/
-		TwoDMatrix&mx = *(new TwoDMatrix(nMax_lag+1, nExog));
-		mx.zeros(); // initialise shocks to 0s
-
-		eval(y, mx, modParams, nMax_lag, residual, g1, g2);	
-};
-
+/*
+* Copyright (C) 2008-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 <http://www.gnu.org/licenses/>.
+*/
+
+/****************************************************** 
+// k_order_perturbation.cpp : Defines the entry point for the k-order perturbation application DLL.
+//
+// called from Dynare dr1_k_order.m, (itself called form resol.m instead of regular dr1.m)
+//            if options_.order < 2 % 1st order only
+//                [ysteady, ghx_u]=k_order_perturbation(dr,task,M_,options_, oo_ , ['.' mexext]);
+//            else % 2nd order
+//                [ysteady, ghx_u, g_2]=k_order_perturbation(dr,task,M_,options_, oo_ , ['.' mexext]);
+// inputs: 
+//			dr,		- Dynare structure
+//			task,  - check or not, not used
+//			M_		- Dynare structure
+//			options_ - Dynare structure
+//			oo_		- Dynare structure 
+//			['.' mexext] Matlab dll extension
+// returns:
+//			 ysteady steady state
+//			ghx_u - first order rules packed in one matrix
+//			g_2 - 2nd order rules packed in one matrix
+**********************************************************/
+
+
+#include "k_ord_dynare.h"
+#include "math.h"
+#include <cstring>
+
+#include <cctype>
+
+#ifdef  _MSC_VER  //&&WINDOWS 
+
+
+BOOL APIENTRY DllMain( HANDLE hModule, 
+                      DWORD  ul_reason_for_call, 
+                      LPVOID lpReserved
+                      )
+{
+  switch (ul_reason_for_call)
+  {
+  case DLL_PROCESS_ATTACH:
+  case DLL_THREAD_ATTACH:
+  case DLL_THREAD_DETACH:
+  case DLL_PROCESS_DETACH:
+    break;
+  }
+  return TRUE;
+}
+
+// Some MS Windows preambles
+// This is an example of an exported variable
+K_ORDER_PERTURBATION_API int nK_order_perturbation=0;
+
+// This is an example of an exported function.
+K_ORDER_PERTURBATION_API int fnK_order_perturbation(void)
+{
+  return 42;
+}
+// This is the constructor of a class that has been exported.
+// see k_order_perturbation.h for the class definition
+CK_order_perturbation::CK_order_perturbation()
+{ 
+  return; 
+}
+
+#endif // _MSC_VER && WINDOWS
+
+extern "C" {
+  
+  // mexFunction: Matlab Inerface point and the main application driver	
+  
+  void mexFunction(int nlhs, mxArray* plhs[],
+    int nrhs, const mxArray* prhs[])
+  {
+    if (nrhs < 5)
+      mexErrMsgTxt("Must have at least 5 input parameters.\n");
+    if (nlhs == 0)
+      mexErrMsgTxt("Must have at least 1 output parameter.\n");
+    
+    const mxArray* dr = prhs[0];
+    const int check_flag = (int)mxGetScalar(prhs[1]);
+    const mxArray* M_ = prhs[2];
+    const mxArray* options_=  (prhs[3]);
+    const mxArray* oo_ =  (prhs[4]);
+    
+    mxArray* mFname =mxGetField(M_, 0, "fname");
+    if(!mxIsChar(mFname)){
+      mexErrMsgTxt("Input must be of type char.");
+    }
+    const char*  fName = mxArrayToString(mFname);
+    const char * dfExt=NULL ;//Dyanamic file extension, e.g.".dll" or .mexw32; 
+    if (prhs[5] != NULL){
+      const mxArray* mexExt = prhs[5];
+      dfExt= mxArrayToString(mexExt);
+    } 
+    
+#ifdef DEBUG		
+    mexPrintf("k_order_perturbation: check_flag = %d ,  fName = %s , mexExt=%s.\n", check_flag,fName, dfExt);
+#endif		
+    int kOrder;
+    mxArray* mxFldp = mxGetField(options_, 0,"order" );
+    if (mxIsNumeric(mxFldp))
+      kOrder = (int)mxGetScalar(mxFldp);
+    else
+      kOrder = 1;
+    
+    double qz_criterium = 1+1e-6;
+    mxFldp = mxGetField(options_, 0,"qz_criterium" );
+    if (mxIsNumeric(mxFldp))
+      qz_criterium = (double)mxGetScalar(mxFldp);
+    
+    mxFldp 	= mxGetField(M_, 0,"params" );
+    double * dparams = (double *) mxGetData(mxFldp);
+    int npar = (int)mxGetM(mxFldp);
+    Vector * modParams =  new Vector(dparams, npar);
+#ifdef DEBUG		
+    mexPrintf("k_ord_perturbation: qz_criterium=%g, nParams=%d .\n",qz_criterium,npar);  
+    for (int i = 0; i < npar; i++) {
+      mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]);  }
+    //	for (int i = 0; i < npar; i++) {
+    //        mexPrintf("k_ord_perturbation: params_vec[%d]= %g.\n", i, params_vec[i] );   }
+#endif		        
+    
+    mxFldp 	= mxGetField(M_, 0,"Sigma_e" );
+    dparams = (double *) mxGetData(mxFldp);
+    npar = (int)mxGetN(mxFldp);
+    TwoDMatrix * vCov =  new TwoDMatrix(npar, npar, dparams);
+    
+    
+    //		mxFldp 	= mxGetField(oo_, 0,"steady_state" ); // use in order of declaration
+    mxFldp 	= mxGetField(dr, 0,"ys" );  // and not in order of dr.order_var
+    //		mxFldp 	= mxGetField(oo_, 0,"dyn_ys" );  // and NOT extended ys
+    dparams = (double *) mxGetData(mxFldp);
+    const int nSteady = (int)mxGetM(mxFldp);
+    Vector * ySteady =  new Vector(dparams, nSteady);
+    
+    
+    mxFldp = mxGetField(dr, 0,"nstatic" );
+    const int nStat = (int)mxGetScalar(mxFldp);
+    mxFldp = mxGetField(dr, 0,"npred" );
+    int nPred = (int)mxGetScalar(mxFldp);
+    mxFldp = mxGetField(dr, 0,"nspred" );
+    const int nsPred = (int)mxGetScalar(mxFldp);
+    mxFldp = mxGetField(dr, 0,"nboth" );
+    const int nBoth = (int)mxGetScalar(mxFldp);
+    mxFldp = mxGetField(dr, 0,"nfwrd" );
+    const int nForw = (int)mxGetScalar(mxFldp);
+    mxFldp = mxGetField(dr, 0,"nsfwrd" );
+    const int nsForw = (int)mxGetScalar(mxFldp);
+    
+    mxFldp = mxGetField(M_, 0,"exo_nbr" );
+    const int nExog = (int)mxGetScalar(mxFldp);
+    mxFldp = mxGetField(M_, 0,"endo_nbr" );
+    const int nEndo = (int)mxGetScalar(mxFldp);
+    mxFldp = mxGetField(M_, 0,"param_nbr" );
+    const int nPar = (int)mxGetScalar(mxFldp);
+    // it_ should be set to M_.maximum_lag
+    mxFldp = mxGetField(M_, 0,"maximum_lag" );
+    const int nMax_lag = (int)mxGetScalar(mxFldp);
+    
+    nPred -= nBoth; // correct nPred for nBoth.
+    
+    mxFldp 	= mxGetField(dr, 0,"order_var" );
+    dparams = (double *) mxGetData(mxFldp);
+    npar = (int)mxGetM(mxFldp);
+    if (npar != nEndo) {  //(nPar != npar)
+      mexErrMsgTxt("Incorrect number of input var_order vars.\n");
+      //return;
+    } 
+    vector<int> * var_order_vp = (new vector<int>(nEndo));//nEndo));
+    for (int v=0;v<nEndo;v++){
+      (*var_order_vp)[v] =(int)(*(dparams++));//[v];
+#ifdef DEBUG		
+      mexPrintf("var_order_vp)[%d] = %d .\n", v, (*var_order_vp)[v] );
+#endif
+    }
+    
+    // the lag, current and lead blocks of the jacobian respectively
+    mxFldp 	= mxGetField(M_, 0,"lead_lag_incidence" );
+    dparams = (double *) mxGetData(mxFldp);
+    npar = (int)mxGetN(mxFldp);
+    int nrows = (int)mxGetM(mxFldp);
+#ifdef DEBUG		
+    mexPrintf("ll_Incidence nrows=%d, ncols = %d .\n", nrows, npar);
+#endif
+    TwoDMatrix * llincidence =  new TwoDMatrix ( nrows, npar,dparams);
+    if (npar != nEndo){
+      mexPrintf("Incorrect length of lead lag incidences: ncol=%d != nEndo =%d .\n", npar, nEndo);
+      return;
+    }
+#ifdef DEBUG		
+    for (int j=0;j<nrows;j++){
+      for (int i=0;i<nEndo;i++){
+        mexPrintf("llincidence->get(%d,%d) =%d .\n", 
+          j, i, (int)llincidence->get(j,i));}}
+#endif
+    
+    const int jcols = nExog+nEndo+nsPred+nsForw; // Num of Jacobian columns
+    mexPrintf("k_order_perturbation: jcols= %d .\n", jcols);
+    
+    mxFldp= mxGetField(M_, 0,"var_order_endo_names" );
+    mexPrintf("k_order_perturbation: Get nendo .\n");
+    const int nendo = (int)mxGetM(mxFldp);
+    const int widthEndo = (int)mxGetN(mxFldp);
+    const char**  endoNamesMX= DynareMxArrayToString( mxFldp,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 = (int)mxGetM(mxFldp);
+    const int widthExog = (int)mxGetN(mxFldp);
+    const char**  exoNamesMX= DynareMxArrayToString( mxFldp,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)
+      mexErrMsgTxt("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[%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 = 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: Calling dynamicDLL constructor.\n");
+#endif				
+      //			DynamicFn * pDynamicFn = loadModelDynamicDLL (fname);
+      DynamicModelDLL dynamicDLL(fName,nEndo, jcols, nMax_lag, nExog, dfExt);
+      
+      // intiate tensor library
+#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 dynare constructor .\n");
+#endif			
+      // make KordpDynare object
+      KordpDynare dynare(endoNamesMX,  nEndo, exoNamesMX,  nExog, nPar, // paramNames,
+        ySteady, vCov, modParams, nStat, nPred, nForw, nBoth,
+        jcols, nSteps, kOrder, journal, dynamicDLL, sstol, var_order_vp, 
+        llincidence,qz_criterium );
+      
+      // construct main K-order approximation class
+#ifdef DEBUG		
+      mexPrintf("k_order_perturbation: Call Approximation constructor with qz_criterium=%f \n", qz_criterium);
+#endif
+      Approximation app(dynare, journal,  nSteps, false, qz_criterium);
+      // run stochastic steady 
+#ifdef DEBUG		
+      mexPrintf("k_order_perturbation: Calling walkStochSteady.\n");
+#endif			
+      app.walkStochSteady();		
+      
+      
+      ConstTwoDMatrix ss(app.getSS());
+      // 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);
+      }
+      
+      std::string ss_matrix_name(fName);//params.prefix);
+      ss_matrix_name += "_steady_states";
+      //			ConstTwoDMatrix(app.getSS()).writeMat4(matfd, ss_matrix_name.c_str());
+      ss.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);
+      
+      /* Write derivative outputs into memory map */
+      map<string,ConstTwoDMatrix> mm;
+      app.getFoldDecisionRule().writeMMap(&mm);
+      
+#ifdef DEBUG		
+      app.getFoldDecisionRule().print();
+      mexPrintf("k_order_perturbation: Map print: \n");
+      for (map<string,ConstTwoDMatrix>::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());
+      
+      
+      // developement of the output.
+#ifdef DEBUG		
+      mexPrintf("k_order_perturbation: Filling 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. */
+        TwoDMatrix tmp_ss(nEndo,1, mxGetPr(plhs[0]));
+        tmp_ss = (const TwoDMatrix&)ss;
+#ifdef DEBUG		
+        //				tmp_ss.print();  !! This print Crashes???
+#endif					
+      }
+      if (nlhs >= 2)
+      {
+        /* Set the output pointer to the combined output matrix gyu = [gy gu]. */
+        int ii=1;
+        for (map<string,ConstTwoDMatrix>::const_iterator cit=mm.begin();
+        ((cit !=mm.end())&&(ii<nlhs)); ++cit) 
+        {
+          //if ((*cit).first!="g_0") 
+          {
+            plhs[ii] = mxCreateDoubleMatrix((*cit).second.numRows(), (*cit).second.numCols(), mxREAL);
+            TwoDMatrix dgyu((*cit).second.numRows(), (*cit).second.numCols(), mxGetPr(plhs[ii]));
+            dgyu = (const TwoDMatrix&)(*cit).second;
+#ifdef DEBUG		
+            mexPrintf("k_order_perturbation: cit %d print: \n", ii);
+            (*cit).second.print();
+            mexPrintf("k_order_perturbation: dguy %d print: \n", ii);
+            //                      dgyu.print(); !! This print Crashes???
+#endif					
+            ++ii;
+          }
+        }
+      }
+    } catch (const KordException& e) {
+      printf("Caugth Kord exception: ");
+      e.print();
+      mexPrintf("Caugth Kord exception: %s", e.get_message());
+      std::string errfile(fName);//(params.basename);
+      errfile += "_error.log";
+      FILE* errfd = NULL;
+      if (NULL == (errfd=fopen(errfile.c_str(), "wb"))) {
+        fprintf(stderr, "Couldn't open %s for writing.\n", errfile.c_str());
+        return;// e.code();
+      }
+      fprintf(errfd ,"Caugth Kord exception: %s", e.get_message());
+      fclose(errfd);
+      return;// e.code();
+    } catch (const TLException& e) {
+      printf("Caugth TL exception: ");
+      e.print();
+      return;// 255;
+    } catch (SylvException& e) {
+      printf("Caught Sylv exception: ");
+      e.printMessage();
+      return;// 255;
+    } catch (const DynareException& e) {
+      printf("Caught KordpDynare exception: %s\n", e.message());
+      mexPrintf("Caugth Dynare exception: %s", e.message());
+      std::string errfile(fName);//(params.basename);
+      errfile += "_error.log";
+      FILE* errfd = NULL;
+      if (NULL == (errfd=fopen(errfile.c_str(), "wb"))) {
+        fprintf(stderr, "Couldn't open %s for writing.\n", errfile.c_str());
+        return;// e.code();
+      }
+      fprintf(errfd ,"Caugth KordDynare  exception: %s", e.message());
+      fclose(errfd);
+      return;// 255;
+    } catch (const ogu::Exception& e) {
+      printf("Caught ogu::Exception: ");
+      e.print();
+      mexPrintf("Caugth general exception: %s", e.message());
+      std::string errfile(fName);//(params.basename);
+      errfile += "_error.log";
+      FILE* errfd = NULL;
+      if (NULL == (errfd=fopen(errfile.c_str(), "wb"))) {
+        fprintf(stderr, "Couldn't open %s for writing.\n", errfile.c_str());
+        return;// e.code();
+      }
+      e.print(errfd);
+      fclose(errfd);
+      return;// 255;
+    }  //catch
+  }; // end of mexFunction()
+}; // end of extern C
+
+//////////////////////////////////////////////////////
+// Convert Matlab Dynare endo and exo names array to C type array of string pointers
+// Poblem is that Matlab mx function returns a long string concatenated by columns rather than rows
+// hence a rather low level approach is needed
+///////////////////////////////////////////////////////
+const char ** DynareMxArrayToString(const mxArray * mxFldp,const int len,const int width )
+{
+  char * cNamesCharStr= mxArrayToString(mxFldp);
+  const char ** ret = DynareMxArrayToString(cNamesCharStr,len, width );
+  return ret;
+}
+
+const char ** DynareMxArrayToString(const char * cNamesCharStr,const int len,const int width )
+{
+  char cNamesMX[len][width+1] ;//
+#ifdef DEBUG
+  mexPrintf("loop DynareMxArrayToString cNamesCharStr = %s \n", cNamesCharStr);
+#endif
+  for (int i=0;i<width;i++){
+    for (int j=0;j<len;j++){
+      // Allow alphanumeric and underscores "_" only:
+      if (isalnum(cNamesCharStr[j+i*len])||('_'==cNamesCharStr[j+i*len])){
+        cNamesMX[j][i]=cNamesCharStr[j+i*len];
+      }
+      else cNamesMX[j][i]='\0';
+    }
+  }
+  const char ** ret= (const char **)mxCalloc (len, sizeof(char*));
+  for (int j=0;j<len;j++){
+				cNamesMX[j][width]='\0';
+#ifdef DEBUG
+        //				mexPrintf("String [%d]= %s \n", j, cNamesMX[j]);
+#endif
+        char * token= (char*) mxCalloc ( strlen(cNamesMX[j])+1,sizeof(char));
+        strcpy (token, cNamesMX[j]);
+        ret[j]=token;
+#ifdef DEBUG
+        mexPrintf("ret [%d]= %s \n", j, ret[j]);
+#endif
+  }
+  return ret;
+}
+
+
+/***********************************
+* Members of DynamicModelDLL for handling loading and calling 
+* <model>_dynamic () function
+**************************************/
+DynamicModelDLL::DynamicModelDLL(const char * modName, const int y_length, const int j_cols, 
+                                 const int n_max_lag, const int n_exog, const char * sExt)
+                                 : length(y_length),jcols( j_cols), nMax_lag(n_max_lag), nExog(n_exog)
+{
+  char fName[MAX_MODEL_NAME];
+  strcpy(fName,modName);
+  using namespace std;
+  strcat(fName,"_dynamic");
+#ifdef DEBUG
+		mexPrintf("MexPrintf: Call Load run DLL %s .\n", fName);
+#endif	
+    
+    try {
+#ifdef WINDOWS
+      if (sExt==NULL) sExt=(".dll");
+      HINSTANCE dynamicHinstance;
+      //		dynamicHinstance=::LoadLibraryEx(strcat(fNname,"_.dll"),NULL,DONT_RESOLVE_DLL_REFERENCES);//sExt); //"_.dll");
+      dynamicHinstance=::LoadLibrary(strcat(fName,sExt));//.dll); //"_.dll");
+      if (dynamicHinstance==NULL)
+        throw 1; //alt: return;
+      //		(DynamicFn*)	typedef void * (__stdcall *DynamicFn)();
+#ifdef DEBUG
+      mexPrintf("MexPrintf: Call GetProcAddress  %s .\n", fName);
+#endif
+      Dynamic = (DynamicFn *) ::GetProcAddress(dynamicHinstance,"Dynamic");
+      
+# else // __linux__
+      if (sExt==NULL) sExt=(".so");
+      dynamicHinstance = dlopen(strcat(fName,sExt), RTLD_NOW);
+      if((dynamicHinstance == NULL) || dlerror()){
+        cerr << dlerror() << endl;
+        mexPrintf("MexPrintf:Error loading DLL: %s", dlerror);
+        throw 1;
+      }
+      Dynamic = (DynamicFn)dlsym(dynamicHinstance, "Dynamic");
+      if((Dynamic  == NULL) || dlerror()){
+        cerr << dlerror() << endl;
+        mexPrintf("MexPrintf:Error finding DLL function: %s", dlerror);
+        throw 2;
+      }
+# endif
+      
+    } catch (int i) {
+      mexPrintf("MexPrintf: error in Load and run DLL %s , %d.\n", fName, i);
+      mexErrMsgTxt("Err: An error in Load and run DLL  .\n");
+      return;
+      
+    } catch (...) {
+      mexPrintf("MexPrintf: Unknown error in Call MATLAB function %s.\n", fName);
+      mexErrMsgTxt("Err: Unknown error in Load and run DLL  .\n");
+      return;
+    }
+}
+
+// close DLL: If the referenced object was successfully closed, 
+// close() returns 0, non 0 otherwise
+int DynamicModelDLL::close(){
+#ifdef WINDOWS
+  // MS FreeLibrary returns non 0 if OK, 0 if fails.
+  bool rb=FreeLibrary(dynamicHinstance);
+  if (rb)
+    return 0;
+  else 
+    return 1;
+# else // linux
+  //If OK, dlclose() returns 0, non 0 otherwise
+  return dlclose(dynamicHinstance);
+# endif
+};
+
+
+void DynamicModelDLL::eval(const Vector&y, const TwoDMatrix&x, const  Vector* modParams, 
+                           int it_, Vector&residual, TwoDMatrix*g1, TwoDMatrix*g2){
+  
+  double  *dresidual, *dg1=NULL, *dg2=NULL; 
+  //int length=y.length(); // not!
+		if ((jcols-nExog)!=y.length()){
+      // throw DLL Error
+      mexPrintf(" DLL Error: (jcols-nExog)!=ys.length() \n");
+      return;
+    }
+    if (residual.length()<length){ // dummy or insufficient
+      Vector*tempv= new Vector(length );
+      residual=*tempv;
+      delete tempv;
+      residual.zeros();
+    }
+    if (g1!=NULL){
+      if (g1->nrows()!=length){ // dummy
+        delete g1;
+        g1=	new TwoDMatrix( length, jcols); // and get a new one
+        g1->zeros();
+      }
+      dg1= const_cast<double*>(g1->base());
+    }
+    if (g2!=NULL){
+      if (g2->nrows()!=length){ // dummy 
+        delete g2;
+        g2=	new TwoDMatrix( length, jcols*jcols);// and get a new one
+        g2->zeros();
+      }
+      dg2= const_cast<double*>(g2->base());
+    }
+    dresidual=const_cast<double*>(residual.base());
+    double *dy=const_cast<double*>(y.base());
+    double *dx=const_cast<double*>(x.base());
+    double *dbParams=const_cast<double*>(modParams->base());
+#ifdef DEBUG		
+    mexPrintf(" try eval Dynamic with ne g1: cols=%d , rows=%d\n"
+      , g1->ncols(),g1->nrows());
+    for (int i = 0; i < modParams->length(); i++) {
+      mexPrintf("k_ord_perturbation: Params[%d]= %g.\n", i, (*modParams)[i]);  }
+    for (int i = 0; i < jcols-nExog; i++) {
+      mexPrintf("k_ord_perturbation: Ys[%d]= %g.\n", i, dy[i]);} 
+    mexPrintf("k_order_perturbation: call <model> Dynamic dParams= %g ,  , dy = %g dx = %f .\n"
+      ,dbParams[0],dy[0],dx[0]);
+    
+#endif        
+    try{
+      Dynamic(dy, dx, nExog, dbParams, it_, dresidual, dg1, dg2);
+    }catch (...){
+      mexPrintf("MexPrintf: error in run Dynamic DLL \n");
+    }
+};
+
+void DynamicModelDLL::eval(const Vector&y, const TwoDMatrix&x, const Vector * modParams, 
+                           Vector&residual, TwoDMatrix*g1, TwoDMatrix*g2){
+  
+		eval(y, x, modParams, nMax_lag, residual, g1, g2);	
+};
+
+void DynamicModelDLL::eval(const Vector&y, const Vector&x, const Vector * modParams, 
+                           Vector&residual, TwoDMatrix*g1, TwoDMatrix*g2){
+  
+                           /** ignore given exogens and create new 2D x matrix since 
+                           * when calling <model>_dynamic(z,x,params,it_) x must be equal to
+                           * zeros(M_.maximum_lag+1,M_.exo_nbr)
+		**/
+		TwoDMatrix&mx = *(new TwoDMatrix(nMax_lag+1, nExog));
+    mx.zeros(); // initialise shocks to 0s
+    
+    eval(y, mx, modParams, nMax_lag, residual, g1, g2);	
+};
+
diff --git a/mex/sources/korderpert/src/k_order_perturbation.h b/mex/sources/korderpert/src/k_order_perturbation.h
index 9e8ee3fff..6f9df47a4 100644
--- a/mex/sources/korderpert/src/k_order_perturbation.h
+++ b/mex/sources/korderpert/src/k_order_perturbation.h
@@ -1,110 +1,108 @@
-/*
- * Copyright (C) 2008-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 <http://www.gnu.org/licenses/>.
- */
-
-// The following ifdef block is the standard way of creating macros which make exporting 
-// from a DLL simpler. All files within this DLL are compiled with the K_ORDER_PERTURBATION_EXPORTS
-// symbol defined on the command line. this symbol should not be defined on any project
-// that uses this DLL. This way any other project whose source files include this file see 
-// K_ORDER_PERTURBATION_API functions as being imported from a DLL, wheras this DLL sees symbols
-// defined with this macro as being exported.
-
-#ifdef WINDOWS 
-#ifdef K_ORDER_PERTURBATION_EXPORTS
-#define K_ORDER_PERTURBATION_API __declspec(dllexport)
-#else
-#define K_ORDER_PERTURBATION_API __declspec(dllimport)
-#endif
-
-#include "stdafx.h"
-
-#else
-#include <dlfcn.h> // unix/linux DLL (.so) handling routines
-#endif
-
-#include <string>
-#include "mex.h"
-
-// <model>_Dynamic DLL pointer
-#ifdef WINDOWS
-    typedef void  *(DynamicFn)
-#else // linux
-    typedef void  (*DynamicFn)
-#endif
-(double *y, double *x, int nb_row_x, double *params, 
- int it_, double *residual, double *g1, double *g2);
-
-//DynamicFn Dynamic;
-
-typedef void *(mexFunctionPtr)(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]);
-
-const int MAX_MODEL_NAME=100;
-
-/**
-* creates pointer to Dynamic function inside <model>_dynamic.dll
-* and handles calls to it.
-**/
-class DynamicModelDLL
-{
-private:
-#ifdef WINDOWS
-        DynamicFn  *Dynamic;// pointer to the Dynamic function in DLL
-#else
-        DynamicFn  Dynamic;// pointer to the Dynamic function in DLL
-#endif
-	const int length;  // tot num vars = Num of Jacobian rows
-	const int jcols;  // tot num var t-1, t and t+1 instances + exogs = Num of Jacobian columns
-	const int nMax_lag; // no of lags
-	const int nExog; // no of exogenous
-//	const char * sExt; // dynamic file extension: dll, mexw32, ...
-#ifdef WINDOWS
-	HINSTANCE dynamicHinstance;  // DLL instance pointer in Windows
-# else // linux
-	void * dynamicHinstance ;	// and in Linux
-#endif
-
-
-public:
-	// construct and load Dynamic model DLL 
-	DynamicModelDLL(const char* fname, const int length,const int jcols, 
-		const int nMax_lag, const int nExog, const char * sExt);
-	virtual ~DynamicModelDLL(){close();};
-	//	DynamicFn get(){return DynamicDLLfunc;};
-	//    void 
-	//    ((DynamicFn())*) get(){return Dynamic;};
-	// evaluate Dynamic model DLL
-	void eval(double *y, double *x, int nb_row_x, double *params, 
-		int it_, double *residual, double *g1, double *g2){
-		Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2);
-	};
-	void eval(const Vector&y, const Vector&x,  const Vector* params, 
-		Vector&residual, TwoDMatrix*g1, TwoDMatrix*g2);
-	void eval(const Vector&y, const TwoDMatrix&x,  const Vector* params, 
-		int it_, Vector&residual, TwoDMatrix*g1, TwoDMatrix*g2);
-	void eval(const Vector&y, const TwoDMatrix&x,  const Vector* params, 
-		Vector& residual, TwoDMatrix *g1, TwoDMatrix *g2);
-//	void eval(const Vector&y, const TwoDMatrix&x,  const Vector* params, 
-//		Vector& residual, double *g1, double *g2);
-	// close DLL: If the referenced object was successfully closed, 
-	// close() returns 0, non 0 otherwise
-	int close();
-	
-};
-// convert Matlab endo and exo names array to C type array of strings
-const char ** DynareMxArrayToString(const mxArray * mxFldp, const int len, const int width );
-const char ** DynareMxArrayToString(const char * cArray, const int len, const int width );
+/*
+* Copyright (C) 2008-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 <http://www.gnu.org/licenses/>.
+*/
+
+// The following ifdef block is the standard way of creating macros which make exporting 
+// from a DLL simpler. All files within this DLL are compiled with the K_ORDER_PERTURBATION_EXPORTS
+// symbol defined on the command line. this symbol should not be defined on any project
+// that uses this DLL. This way any other project whose source files include this file see 
+// K_ORDER_PERTURBATION_API functions as being imported from a DLL, wheras this DLL sees symbols
+// defined with this macro as being exported.
+
+#ifdef WINDOWS 
+#ifdef K_ORDER_PERTURBATION_EXPORTS
+#define K_ORDER_PERTURBATION_API __declspec(dllexport)
+#else
+#define K_ORDER_PERTURBATION_API __declspec(dllimport)
+#endif
+
+#include "stdafx.h"
+
+#else
+#include <dlfcn.h> // unix/linux DLL (.so) handling routines
+#endif
+
+#include <string>
+#include "mex.h"
+
+// <model>_Dynamic DLL pointer
+#ifdef WINDOWS
+typedef void  *(DynamicFn)
+#else // linux
+typedef void  (*DynamicFn)
+#endif
+(double *y, double *x, int nb_row_x, double *params, 
+ int it_, double *residual, double *g1, double *g2);
+
+//DynamicFn Dynamic;
+
+typedef void *(mexFunctionPtr)(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]);
+
+const int MAX_MODEL_NAME=100;
+
+/**
+* creates pointer to Dynamic function inside <model>_dynamic.dll
+* and handles calls to it.
+**/
+class DynamicModelDLL
+{
+private:
+#ifdef WINDOWS
+  DynamicFn  *Dynamic;// pointer to the Dynamic function in DLL
+#else
+  DynamicFn  Dynamic;// pointer to the Dynamic function in DLL
+#endif
+  const int length;  // tot num vars = Num of Jacobian rows
+  const int jcols;  // tot num var t-1, t and t+1 instances + exogs = Num of Jacobian columns
+  const int nMax_lag; // no of lags
+  const int nExog; // no of exogenous
+  //	const char * sExt; // dynamic file extension: dll, mexw32, ...
+#ifdef WINDOWS
+  HINSTANCE dynamicHinstance;  // DLL instance pointer in Windows
+# else // linux
+  void * dynamicHinstance ;	// and in Linux
+#endif
+  
+  
+public:
+  // construct and load Dynamic model DLL 
+  DynamicModelDLL(const char* fname, const int length,const int jcols, 
+    const int nMax_lag, const int nExog, const char * sExt);
+  virtual ~DynamicModelDLL(){close();};
+
+  // evaluate Dynamic model DLL
+  void eval(double *y, double *x, int nb_row_x, double *params, 
+    int it_, double *residual, double *g1, double *g2){
+    Dynamic(y, x, nb_row_x, params, it_, residual, g1, g2);
+  };
+  void eval(const Vector&y, const Vector&x,  const Vector* params, 
+    Vector&residual, TwoDMatrix*g1, TwoDMatrix*g2);
+  void eval(const Vector&y, const TwoDMatrix&x,  const Vector* params, 
+    int it_, Vector&residual, TwoDMatrix*g1, TwoDMatrix*g2);
+  void eval(const Vector&y, const TwoDMatrix&x,  const Vector* params, 
+    Vector& residual, TwoDMatrix *g1, TwoDMatrix *g2);
+  //	void eval(const Vector&y, const TwoDMatrix&x,  const Vector* params, 
+  //		Vector& residual, double *g1, double *g2);
+  // close DLL: If the referenced object was successfully closed, 
+  // close() returns 0, non 0 otherwise
+  int close();
+  
+};
+// convert Matlab endo and exo names array to C type array of strings
+const char ** DynareMxArrayToString(const mxArray * mxFldp, const int len, const int width );
+const char ** DynareMxArrayToString(const char * cArray, const int len, const int width );
diff --git a/mex/sources/korderpert/src/k_order_test_main.cpp b/mex/sources/korderpert/src/k_order_test_main.cpp
index faf2509f5..8c2311755 100644
--- a/mex/sources/korderpert/src/k_order_test_main.cpp
+++ b/mex/sources/korderpert/src/k_order_test_main.cpp
@@ -1,324 +1,324 @@
-/*
- * Copyright (C) 2008-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 <http://www.gnu.org/licenses/>.
- */
-
-/*************************************
-* 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 "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=NULL ;//Dyanamic 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);
-
-
-		//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<int> * var_order_vp = new vector<int>(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)
-			mexErrMsgTxt("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, 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<string,ConstTwoDMatrix> mm;
-			app.getFoldDecisionRule().writeMMap(&mm);
-#ifdef DEBUG		
-			app.getFoldDecisionRule().print();
-			mexPrintf("k_order_perturbation: Map print: \n");
-			for (map<string,ConstTwoDMatrix>::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;
-}
+/*
+* Copyright (C) 2008-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 <http://www.gnu.org/licenses/>.
+*/
+
+/*************************************
+* 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 "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=NULL ;//Dyanamic 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);
+    
+    
+    //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<int> * var_order_vp = new vector<int>(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)
+      mexErrMsgTxt("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, 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<string,ConstTwoDMatrix> mm;
+      app.getFoldDecisionRule().writeMMap(&mm);
+#ifdef DEBUG		
+      app.getFoldDecisionRule().print();
+      mexPrintf("k_order_perturbation: Map print: \n");
+      for (map<string,ConstTwoDMatrix>::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;
+}