function dW = findxie(x, D, Y, N, K, C, L, Lp, options, w, z, p) 

global oo_

p.xie = x; 

Xnew         = fsolve('findequilibrium', [D; Y; N], options, w, z, p, 'market', 'new');

[~, qnew, ~, ~, Znew, Wnew, Lpnew, Knew, Bnew, Unew, ~, Cnew, Munew] = findequilibrium(Xnew, w, z, p, 'market', 'new');


% New Steady State Allocations

Dnew          = Xnew(1); 
Ynew          = Xnew(2); 
Nnew          = Xnew(3); 


Gnew          = (1 + p.xie)*(1 + p.xis - 1/Munew);
Qnew          = 1/(1 - p.beta*(1 - p.varphi))*Gnew*Ynew/Nnew; 

Lnew          = Lpnew + p.F*p.varphi*Nnew; 
Rnew          = 1/p.beta - 1 + p.delta; 

ynew          = qnew*Ynew; 




% Now start doing transition dynamics

% See if we can use dynare: parameterize Z(N) with chebyshev polynomials

global cfun fspace; 

npts      = 11; 

optset('chebnode', 'nodetype', 2);              % Lobatto nodes (include endpoints)

fspace    = fundefn('cheb', npts, 1, 2); 

Ngrid     = funnode(fspace); 


Zgrid     = zeros(npts, 1); 
Mugrid    = zeros(npts, 1); 

options.Display = 'none';

for i = 1 : npts

    D      = fsolve('findequilibrium', D, options, w, z, p, 'market', 'old', Ngrid(i));

    [~, ~, ~, ~, Zgrid(i), ~, ~, ~, ~, ~, ~, ~, Mugrid(i)]  = findequilibrium(D, w, z, p, 'market', 'old', Ngrid(i));
    
end


Ggrid  = (1 + p.xie)*(1 + p.xis - 1./Mugrid);               % flow gains from varieties

cfun   = funfitxy(fspace, Ngrid, [Zgrid, Ggrid, Mugrid]);

xpar   = [p.beta; p.varphi; p.delta; p.psi; p.F; p.nu; p.alpha; p.theta; p.phi];                              save xpar xpar; 

x1_ss  = [K;    N];                                                                                           save x1_ss x1_ss; 
x2_ss  = [Ynew; Cnew; Bnew; Lpnew; Knew; Nnew; Wnew; Rnew; Qnew; Znew; Gnew; Znew/Munew; Munew];              save x2_ss x2_ss; 


%dynare transition noclearall;
dynare transition noclearall nopreprocessoroutput;

% Cleanup

filename = 'transition';

eval(['delete ' filename '*.m']);
eval(['delete ' filename '*.mat']);
eval(['delete ' filename '*.swp']);
eval(['delete ' filename '*.log']);
eval(['delete ' filename '*.eps']);
eval(['delete ' filename '*.pdf']);
eval(['delete ' filename '*.fig']);
eval(['delete ' filename '*.jnl']);

Ct     = oo_.endo_simul(2,  :)';
Lpt    = oo_.endo_simul(4,  :)';
Kt     = oo_.endo_simul(5,  :)';
Nt     = oo_.endo_simul(6,  :)';


% Populate first period values from original steady state

Ct(1)  = C; 
Lpt(1) = Lp; 
Kt(1)  = K; 
Nt(1)  = N; 


Lt     = [L; Lpt(2:end) + p.F*(Nt(2:end) - (1 - p.varphi)*Nt(1 : end-1))]; 


T    = length(Ct) - 1; 


Wold = 1/(1 - p.beta)*(log(C)  - p.psi/(1 + p.nu)*L^(1 + p.nu)); 

Wnew = (p.beta.^(0 : 1 : T - 1))*(log(Ct(2 : end)) - p.psi/(1 + p.nu)*Lt(2:end).^(1 + p.nu)) + p.beta^T/(1 - p.beta)*(log(Ct(end)') - p.psi/(1 + p.nu)*Lt(end)'.^(1 + p.nu));

dW   = (exp((1 - p.beta)*(Wnew - Wold)) - 1)*100;


fprintf('\n')
fprintf('Subsidy, Welfare Gains                       = %9.5f %9.5f \n',  [p.xie, dW]);
fprintf('\n')