diff --git a/mex/sources/korderpert/src/first_order.m b/mex/sources/korderpert/src/first_order.m
new file mode 100644
index 000000000..20e7e1414
--- /dev/null
+++ b/mex/sources/korderpert/src/first_order.m
@@ -0,0 +1,110 @@
+function [gy]=first_order(M_, dr, jacobia)
+% fd = jacobia_
+% reorder jacobia_
+[fd]=k_reOrderedJacobia(M_, jacobia)
+
+%ypart=dr;
+ypart.ny=M_.endo_nbr;
+ypart.nyss=dr.nfwrd+dr.nboth;
+ypart.nys=dr.npred;
+ypart.npred=dr.npred-dr.nboth;
+ypart.nboth=dr.nboth;
+ypart.nforw=dr.nfwrd;
+ypart.nstat =dr.nstatic
+nu=M_.exo_nbr
+
+	off= 1; % = 0 in C
+	fyplus = fd(:,off:off+ypart.nyss-1);
+	off= off+ypart.nyss;
+	fyszero= fd(:,off:off+ypart.nstat-1);
+	off= off+ypart.nstat;
+    fypzero= fd(:,off:off+ypart.npred-1);
+	off= off+ypart.npred;
+	fybzero= fd(:,off:off+ypart.nboth-1);
+	off= off+ypart.nboth;
+	fyfzero= fd(:,off:off+ypart.nforw-1);
+	off= off+ypart.nforw;
+	fymins= fd(:,off:off+ypart.nys-1);
+	off= off+ypart.nys;
+	fuzero= fd(:,off:off+nu-1);
+	off=off+ nu;
+	
+	n= ypart.ny+ypart.nboth;
+	%TwoDMatrix 
+    matD=zeros(n,n);
+%	matD.place(fypzero,0,0);
+	matD(1:n-ypart.nboth,1:ypart.npred)= fypzero;
+%	matD.place(fybzero,0,ypart.npred);
+	matD(1:n-ypart.nboth,ypart.npred+1:ypart.npred+ypart.nboth)=fybzero;
+%	matD.place(fyplus,0,ypart.nys()+ypart.nstat);
+	matD(1:n-ypart.nboth,ypart.nys+ypart.nstat+1:ypart.nys+ypart.nstat+ypart.nyss)=fyplus;
+	for i=1:ypart.nboth
+		matD(ypart.ny()+i,ypart.npred+i)= 1.0;
+    end
+
+    matE=[fymins, fyszero, zeros(n-ypart.nboth,ypart.nboth), fyfzero; zeros(ypart.nboth,n)];
+%	matE.place(fymins;
+%	matE.place(fyszero,0,ypart.nys());
+%	matE.place(fyfzero,0,ypart.nys()+ypart.nstat+ypart.nboth);
+    
+	for i= 1:ypart.nboth
+		matE(ypart.ny()+i,ypart.nys()+ypart.nstat+i)= -1.0;
+    end
+    matE=-matE; %matE.mult(-1.0);
+	
+%    vsl=zeros(n,n);
+%	vsr=zeros(n,n);
+%	lwork= 100*n+16;
+%	work=zeros(1,lwork);
+%	bwork=zeros(1,n);
+	%int info;
+
+%    	LAPACK_dgges("N","V","S",order_eigs,&n,matE.getData().base(),&n,
+%		matD.getData().base(),&n,&sdim,alphar.base(),alphai.base(),
+%		beta.base(),vsl.getData().base(),&n,vsr.getData().base(),&n,
+%		work.base(),&lwork,&(bwork[0]),&info);
+
+   [matE1,matD1,vsr,sdim,dr.eigval,info] = mjdgges(matE,matD,1);
+
+	bk_cond= (sdim==ypart.nys);
+	
+%  	ConstGeneralMatrix z11(vsr,0,0,ypart.nys(),ypart.nys());
+	 z11=vsr(1:ypart.nys,1:ypart.nys);
+%	ConstGeneralMatrix z12(vsr,0,ypart.nys(),ypart.nys(),n-ypart.nys());
+	 z12=vsr(1:ypart.nys(),ypart.nys+1:end);%, n-ypart.nys);
+%	ConstGeneralMatrix z21(vsr,ypart.nys(),0,n-ypart.nys(),ypart.nys());
+	 z21=vsr(ypart.nys+1:end,1:ypart.nys);
+%	ConstGeneralMatrix z22(vsr,ypart.nys(),ypart.nys(),n-ypart.nys(),n-ypart.nys());
+	 z22=vsr(ypart.nys+1:end,ypart.nys+1:end);
+
+% 	GeneralMatrix sfder(z12,"transpose");
+    sfder=z12';%,"transpose");
+%	z22.multInvLeftTrans(sfder);
+    sfder=z22'\sfder;
+	sfder=-sfder;% .mult(-1);
+
+    %s11(matE,0,0,ypart.nys(),ypart.nys());
+     s11=matE1(1:ypart.nys,1:ypart.nys);
+%	 t11=(matD1,0,0,ypart.nys(),ypart.nys());
+	 t11=matD1(1:ypart.nys,1:ypart.nys);
+	 dumm=(s11');%,"transpose");
+	%z11.multInvLeftTrans(dumm);
+    dumm=z11'\dumm;
+	 preder=(dumm');%,"transpose");
+	%t11.multInvLeft(preder);
+    preder=t11\preder;
+	%preder.multLeft(z11);
+    preder= z11*preder;
+	
+%	gy.place(preder,ypart.nstat,0);
+%	gy=(zeros(ypart.nstat,size(preder,2)) ;preder);
+%	 sder(sfder,0,0,ypart.nstat,ypart.nys());
+	 sder=sfder(1:ypart.nstat,1:ypart.nys);
+%	gy.place(sder,0,0);
+%	gy(1:ypart.nstat, 1:ypart.nys)=sder;
+%    gy=[sder;preder];
+%	 fder(sfder,ypart.nstat+ypart.nboth,0,ypart.nforw,ypart.nys());
+	 fder=sfder(ypart.nstat+ypart.nboth+1:ypart.nstat+ypart.nboth+ypart.nforw,1:ypart.nys);
+%	gy.place(fder,ypart.nstat+ypart.nys(),0);
+%	gy(ypart.nstat+ypart.nys,1)=fder;
+    gy=[sder;preder;fder];
diff --git a/mex/sources/korderpert/src/fs2000k.mod b/mex/sources/korderpert/src/fs2000k.mod
index 9ba2b6061..731a946c7 100644
--- a/mex/sources/korderpert/src/fs2000k.mod
+++ b/mex/sources/korderpert/src/fs2000k.mod
@@ -17,9 +17,13 @@
 // Note that there is an initial minus sign missing in equation (A1), p. S63.
 //
 // Michel Juillard, February 2004
+// Modified for testing k_order_perturbation by GP, Jan-Feb 09
+
 options_.usePartInfo=0;
-//var m P c e W R k d n l gy_obs gp_obs Y_obs P_obs y dA P2 c2;
-var m P c e W R k d n l gy_obs gp_obs y dA P2 c2;
+options_.use_k_order=0;
+
+//var m m_1 P P_1 c e W R k d n l gy_obs gp_obs Y_obs P_obs y dA P2 c2;
+var m m_1 P P_1 c e W R k d n l gy_obs gp_obs y dA P2 c2;
 varexo e_a e_m;
 
 parameters alp bet gam mst rho psi del;
@@ -34,7 +38,7 @@ del = 0.02;
 
 model (use_dll);
 dA = exp(gam+e_a);
-log(m) = (1-rho)*log(mst) + rho*log(m(-1))+e_m;
+log(m) = (1-rho)*log(mst) + rho*log(m_1(-1))+e_m;
 -P/(c(+1)*P(+1)*m)+bet*P(+1)*(alp*exp(-alp*(gam+log(e(+1))))*k^(alp-1)*n(+1)^(1-alp)+(1-del)*exp(-(gam+log(e(+1)))))/(c2(+1)*P2(+1)*m(+1))=0;
 W = l/n;
 -(psi/(1-psi))*(c*P/(1-n))+l/n = 0;
@@ -46,16 +50,20 @@ m-1+d = l;
 e = exp(e_a);
 y = k(-1)^alp*n^(1-alp)*exp(-alp*(gam+e_a));
 gy_obs = dA*y/y(-1);
-gp_obs = (P/P(-1))*m(-1)/dA;
+gp_obs = (P/P_1(-1))*m_1(-1)/dA;
 //Y_obs/Y_obs(-1) = gy_obs;
 //P_obs/P_obs(-1) = gp_obs;
 P2 = P(+1);
 c2 = c(+1);
+m_1 = m;
+P_1 = P;
 end;
 
 initval;
 m = mst;
+m_1=mst;
 P = 2.25;
+P_1 = 2.25;
 c = 0.45;
 e = 1;
 W = 4;
@@ -64,11 +72,9 @@ k = 6;
 d = 0.85;
 n = 0.19;
 l = 0.86;
+y = 0.6;
 gy_obs = exp(gam);
 gp_obs = exp(-gam); 
-//  Y_obs = 20000;
-//  P_obs = 51;
-y = 0.6;
 dA = exp(gam);
   P2=P;
   c2=c;
@@ -100,12 +106,13 @@ end;
 //varobs P_obs Y_obs;
 varobs gp_obs gy_obs;
 
+steady(solve_algo = 2);
+
 //observation_trends;
 //P_obs (log(mst)-gam);
 //Y_obs (gam);
 //end;
 
-
 //options_.useAIM = 1;
 estimation(datafile=fsdat,nobs=192,loglinear,mh_replic=2000,
 	mode_compute=4,mh_nblocks=2,mh_drop=0.45,mh_jscale=0.65);