function System06_16_10=System06_16_10(y)

load savedParameters

P=y(1);
lambda=y(2);
M=y(3);


%savedParameters=[epsilon exponent delta deltad alpha chi sigma eta phi beta R rho thetac thetaf gamma A0 kappa0 M0 sig2eps sig2kappa D];

epsilon=savedParameters(1);
exponent=savedParameters(2);
delta=savedParameters(3);
deltad=savedParameters(4);
alpha=savedParameters(5);
chi=savedParameters(6);
sigma=savedParameters(7);
eta=savedParameters(8);
phi=savedParameters(9);
beta=savedParameters(10);
R=savedParameters(11);
rho=savedParameters(12);
thetac=savedParameters(13);
thetaf=savedParameters(14);
gamma=savedParameters(15);
A0=savedParameters(16);
kappa0=savedParameters(17);
M0=savedParameters(18);
sig2eps=savedParameters(19);
sig2kappa=savedParameters(20);
D=savedParameters(21);


xi=kappa0/(1-beta)/M;
C=(P*xi)^(-1/gamma);

N=((R-1+delta)*lambda/epsilon)^(1/exponent);             

Pd=(1-delta)*lambda*sigma*P/(sigma-1);
Pn=((P^(1-sigma)-D*Pd^(1-sigma))/(1-D))^(1/(1-sigma));
Ln=(C/A0*(Pn/P)^(-sigma))^(1/alpha);

w=(sigma-1)*P^(sigma-sigma/alpha)/sigma*alpha*A0^(1/alpha)*C^(1-1/alpha)*Pn^(1+sigma/alpha-sigma);
Ld=(w/P/A0/lambda/alpha)^(1/(alpha-1));
L=D*Ld+(1-D)*Ln;

V=P*C-w*L;
RE=(1-deltad)/deltad*V;
Div=V;



x1=eta*chi*L^phi+(1-eta)*w*xi;
x2=N-(1-delta)/delta*(A0*Ld^alpha-C*(Pd/P)^(-sigma));
x3=M0-M-RE;



x=[x1 x2 x3];

System06_16_10=x;