function [time,Cpath,Kpath,Lpath,etaLpath] = Function15_02_TP(PAR,EXO,VAR)

% Parameter values:
alpha   = PAR.alpha;           
epsilon = PAR.epsilon;
delta   = PAR.delta;    
rho     = PAR.rho; 
mu      = PAR.mu;
Z_0     = PAR.Z_0;



% Exogenous variables:

time      = EXO.time;
t_L       = EXO.t_L;
t_C       = EXO.t_C;
G         = EXO.G;

% Capital
Keq0    = VAR.Keq0;
Keq1    = VAR.Keq1;
Ceq1    = VAR.Ceq1;

% Initial consumption
Ceq0    = VAR.Cbase;


solinit   = bvpinit(linspace(0,400),[Ceq0 Keq0]);
sol       = bvp4c(@capck,@bcfun,solinit);
y         = deval(sol,time);
Cpath     = y(1,:);
Kpath     = y(2,:);


% Dynamics of C and K

function dydx = capck(x,y)

L     = fzero(@(u) eqemp(u,y),0.3);
Y     = Z_0 * L^(1-alpha) * y(2)^alpha ;
dydx(1,1) = (alpha * Y / y(2) - delta - rho) * y(1) ...
    - epsilon * mu * (rho + mu) * y(2) / (1 +t_C) ;
dydx(2,1) = Y - y(1)  - delta * y(2) ;

end


% Boundary conditions

function res = bcfun(ya,yb)
   res = [ya(2) - Keq0
          yb(2) - Keq1];
end


% Equilibrium employment

function z=eqemp(u,x)
   z = (1 - alpha) * (1 - t_L) * Z_0 * x(2)^alpha * (1 - u)  ...
       - ((1-epsilon)/epsilon) * (1 + t_C) * x(1) * u^alpha ;
end


end