function [k,l,alphah,pf,pn,w,im,yh,cf,cn,nu,gn,gf,pg,pro,rm,d,b,ppf,q,ph,lh,nl,chf,po,phi,chi,val1,val2,val3,val4]=myfunction(yn,p,g,c,v,zh,zn,alpha1,alpha2,epsis,mus,deltag,theta,theta2,varphi,tauf,tau,i,a,r,change,rp,rppf,psi,epsilon,yf,sigma,h,delta)
% First group of variables ***********************************************
  % definition of temporary parameters
coef1=alpha2/(alpha1+alpha2);
coef2=-alpha1/(alpha1+alpha2);
coef3=(alpha1/alpha2)^(-coef2);
coef4=(alpha1/alpha2)^coef2;
int1=(alpha1/alpha2)^(alpha2/(alpha1+alpha2))+(alpha1/alpha2)^(-alpha2/(alpha1+alpha2));
int2=(1/(alpha1+alpha2))*int1*(theta2/(theta2-1));
%int2=int1*(theta2/(theta2-1));
  % end of definition of temporary parameters******************************
k=epsis*mus*(1/deltag)*g;

% we solve the non-linear system whom yields labour l and calibrates
% parameter alphah to fit the steady state
G=@(yy)[yn-1+zh*(((1-delta)*yy(1))^(yy(2)))*(k^(1-yy(2)));
    -log(delta)-log(yy(1))+log(coef3)-coef1*( log(yy(2))+log(v)+log(1-yn)-log(1-delta)-log(yy(1)) )+(1/(alpha1+alpha2))*log((yn/(zn*(k^(1-alpha1-alpha2)))));];
options=optimset('Display','Final','TolX',1e-10,'TolFun',1e-10);yy0=[1;0.5];
[yy,val1]=fsolve(G,yy0,options);
l=yy(1);
alphah=yy(2);
lh=(1-delta)*l;
nl=delta*l;
% we solve the non-linear system whom yields prices ph,pn and the wage w
F=@(x)[
    -p+((1-varphi)*(x(1)^(1-theta))+ varphi*(x(2)^(1-theta)))^(1/(1-theta));
    -log(x(3))+log(alphah)+log(v)+log(x(1))+log(1-yn)-log(1-delta)-log(l);
   -log(yn)+log(int2)-coef2*log(delta*l)-coef2*log(x(3))+coef2*log(x(2))-coef1*log(x(2))+coef1*log(x(1))+coef1*( log(coef4)-coef2*log(x(3))+coef2*log(x(1))+(1/(alpha1+alpha2))*log((yn/(zn*(k^(1-alpha1-alpha2))))) )];
x0=[1;1;0.8];
option1=optimset('Display','Final','TolX',1e-15,'TolFun',1e-15);
%option=optimoptions('fsolve','algorithm','Levenberg-Marquardt');
[x,val2]=fsolve(F,x0,option1);
pf=x(1);
pn=x(2);
w=x(3);
im=coef4*(w^(-coef2))*(pf^coef2)*((yn/(zn*(k^(1-alpha1-alpha2))))^(1/(alpha1+alpha2)));
yh=zh*(((1-delta)*l)^(alphah))*(k^(1-alphah));

%we fix the values of consumption ch and cn
cf=(1-varphi)*((pf/p)^(-theta))*c;
cn=varphi*((pn/p)^(-theta))*c;

%Now we calibrate the value of nu to fit the steady state and set the
%corresponded value of gn and gh
gn=yn-cn;
var1=(gn/g)*(pn^theta);
H=@(t)(t*(((1-t)*(pf^(1-theta))+ t*(pn^(1-theta)))^(theta/(1-theta)))-var1);
t0=0.5;
option3=optimset('Display','Final','TolX',1e-15,'TolFun',1e-15);
%options=optimoptions('fsolve','algorithm','Levenberg-Marquardt');
[nu,val3]=fsolve(H,t0,option3);
pg=((1-nu)*(pf^(1-theta))+ nu*(pn^(1-theta)))^(1/(1-theta));
gf=(1-nu)*((pf/pg)^(-theta))*g;
%gv=nu*((pn/pg)^(-theta))*g;

% we define allthing about prices
ppf=pf*(psi/epsilon);
q=epsilon*(ppf/p);
ph=v*pf;
% Now we use budget constraint of agents,governement and foreign reserves
% equation to conclude
pro=v*(pf/p)*(1-yn)+(pn/p)*yn-w*l-(pf/p)*im;
I=@(y)[-change*(1-(1/rppf))+y(2)+y(1)-y(2)*(r/rppf)+a-(cf+im+gf-v*(1-yn))-tauf*pro;
    tau*w*l+q*a+y(2)*(1-(r/rppf))+y(3)*(1-(i/rp))-pg*g;
    (1-tau)*w*l+1.01*(1-tauf)*pro-y(3)*(1-(i/rp))+q*y(1)-c];
option4=optimset('Display','Final','TolX',1e-15,'TolFun',1e-15);
y0=[0;0;0;0];
[y,val4]=fsolve(I,y0,option4);
rm=y(1);
d=y(2);
b=y(3);

%external sector
chf=yh;
phi=(chf/yf)*((ph/(epsilon*ppf))^theta2);
po=ph;

%calibrage de chi *************
chi=(log(w)+log(1-tau)-sigma*log(c)-log(h))/log(l);
%********************************************************
%[k,l,alphah,pf,pn,w,im,yh,cf,cn,nu,gn,gf,pg,pro,rm,d,b,ppf,q,ph,lh,nl,chf,po,phi,chi,val1,val2,val3,val4]=myfunction(yn,p,g,c,v,zh,zn,alpha1,alpha2,epsis,mus,deltag,theta,theta2,varphi,tauf,tau,i,a,r,change,rp,rppf,psi,epsilon,yf,sigma,h,delta)
%[k,l,alphah,pf,pn,w,im,yh,cf,cn,nu,gn,gf,pg,pro,rm,d,b,ppf,q,ph,lh,nl,chf,po,phi,chi,val1,val2,val3,val4]=myfunction(0.65,1,0.21,0.65,1,1,1,0.5,0.45,0.12,0.25,0.1,12,12,0.75,0.3,0.175,1.05,0.15,1.065,0.12,1,1,1,1,1500,2,0.15,0.6)
end