function [err, W, L, Lp, K, B, C, Q, Z, Mu, G]  = findequilibrium(X, p)


Y       = X(1); 
N       = X(2); 
    

Mu      = interpv(N, 1)/(1 + p.xis);        % not adjusted for subsidy, so do it here
Z       = interpv(N, 2); 
G       = interpv(N, 3)*(1 + p.xis); 


R       = p.r + p.delta;

Omega   = Z/Mu; 

W       = (1 - p.alpha)*(((Omega^(1 - p.theta) - (1 - p.phi))/p.phi)^(1/(1 - p.theta))*(R/p.alpha)^(-p.alpha))^(1/(1 - p.alpha));        % Invert the CES Cost aggregator                                              % Labor only factor of production                

Lp      = (1 - p.alpha)*p.phi/Mu*((R/p.alpha)^p.alpha*(W/(1 - p.alpha))^(1 - p.alpha)/Omega)^(1 - p.theta)*Y/W; 
K       =       p.alpha*p.phi/Mu*((R/p.alpha)^p.alpha*(W/(1 - p.alpha))^(1 - p.alpha)/Omega)^(1 - p.theta)*Y/R;
B       = (1 - p.phi)/Mu*(1/Omega)^(1 - p.theta)*Y;

C       = Y - p.delta*K - B;

L       = Lp + p.F*p.varphi*N; 

Q       = Y*G/(1 - p.beta*(1 - p.varphi));

err     = zeros(2, 1); 


err(1)  = W      -   p.psi*C*L^p.nu; 
err(2)  = p.F*W  -  (1 + p.xie)*p.beta*Q;  



