function dW = findxie(x, Y, N, K, C, L, options, p, cfun, fspace) 

global oo_

p.xie        = x; 

Xnew         = fsolve('findequilibrium', [Y; N], options, p);

% Find New Steady State

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

[~, Wnew, Lnew, Lpnew, Knew, Bnew, Cnew, Qnew, Znew, Munew, Gnew]  = findequilibrium(Xnew,   p);

Rnew           = p.r + p.delta; 

% Transition dynamics

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

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


dynare transition noclearall;

% 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,  :)';
Lt     = oo_.endo_simul(4,  :)';

% Populate first period values from original steady state

Ct(1)  = C; 
Lt(1)  = L; 

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')