%--- compute V2
C2 = dC2 + C_ss;
L2 = dL2 + L_ss;
Z2 = dZ2 + Z_ss;
Z2b = [Z_ss2; Z2(1:end-1)];

uu = ufun(C2, L2, Z2b, muM1);
V2 = zeros(N,1);
for i = 1:N
    [ee, V2(i)] = findCCV(0, zeros(N,1), i, ufun, C2, L2, Z2b, bbeta, muM1);
end

%--- compute lambdaC
C1 = dC1 + C_ss1;
L1 = dL1 + L_ss1;
Z1 = dZ1 + Z_ss1;
Z1b = [Z_ss1; Z1(1:end-1)];
lamvec = zeros(irf_length,1);
x0 = lam_ss;
for i = irf_length:-1:1
    lamvec(i) = fsolve(@(x) findCCV(x, V2, i, ufun, C1, L1, Z1b, bbeta, muM1), x0, options);
    x0 = lamvec(i);
end

V1 = zeros(N,1);
for i = 1:irf_length
    [ee, V1(i)] = findCCV(0, zeros(N,1), i, ufun, C1, L1, Z1b, bbeta, muM1);
end