sim1_linear.m seems to return wrong results. The following mod-file calls different solvers in sequence, using the previous result as the respective starting value.

//Endogenous variables
var piV rV yV rrstarV zV Rshock;

//Exogenous variables
varexo epsFG shock epsR;       

//Parameters
parameters alpha beta eta phipi phiy rho sig lambda; // 

//Initialization of parameter values
beta    =0.99;      //Discount factor (=inverse nominal interest rate)
alpha   =0.033;     //sensivity of inflation to output gap, parameter determining marginal cost, 
eta     =1;         //sensivity of inflation to output gap, parameter determining price rigidity
phipi   =1.5;       //parameter determining monetary reaction to inflation
phiy    =0.125;     //parameter determining monetary reaction to output
rho     =0.9;       //AR(1) natural rate
sig     =1;         //intertemporal rate of substitution


lambda=0.25;        // prob. of new info arrival

//model equations
model;

// z Variable
zV = alpha*(yV-yV(-1))+piV;//(pV-pV(-1));

piV= (lambda*alpha/ (1-lambda)) *yV 
+lambda*EXPECTATION(-1)(zV)
+lambda*(1-lambda)*EXPECTATION(-2)(zV)
+lambda*(1-lambda)^2*EXPECTATION(-3)(zV)
+lambda*(1-lambda)^3*EXPECTATION(-4)(zV)
+lambda*(1-lambda)^4*EXPECTATION(-5)(zV)
+lambda*(1-lambda)^5*EXPECTATION(-6)(zV)
+lambda*(1-lambda)^6*EXPECTATION(-7)(zV)
+lambda*(1-lambda)^7*EXPECTATION(-8)(zV)
+lambda*(1-lambda)^8*EXPECTATION(-9)(zV)
;                      //marginal cost is proportional to gap, gap is driver of inflation

//piV=0.99*piV(+1)+0.033*yV;

//IS curve (Output Euler) from Kiley code
yV = yV(+1)-sig*(rV-piV(+1) - rrstarV);

//monetary policy' reaction function
rV = min((1-shock), phipi*(piV(+1) + phiy*yV + epsFG));
//rV = phipi*piV + phiy*yV+Rshock;

Rshock=0.5*Rshock(-1)+epsR;

//Natural real interest rate
rrstarV = rho*rrstarV(-1);
end;

initval;
rV=0;
yV=0;
rrstarV=0;
piV=0;
zV=0;
end;

steady;
check;

shocks;
var shock;
periods 1:15;
values 1;

var epsFG;
periods 16 ;
values -0.0117;
end;

perfect_foresight_setup(periods=30);

perfect_foresight_solver(solve_algo=0); 
nonlinear=oo_.endo_simul;
options_.linear=1;
perfect_foresight_solver(solve_algo=0); 
linear=oo_.endo_simul;
perfect_foresight_solver(solve_algo=0); 
options_.linear=0;
perfect_foresight_solver(solve_algo=0); 
perfect_foresight_solver(solve_algo=4); 

But the result of the nonlinear solvers has non-zero residuals in the linear solver and vice versa, while the nonlinear solvers seem to have consistent results.