function F = solve_rbeqpen(x, alffa, betta, delta, ret, J, mort, n, N_rel_all, N_rel)
%% Initialize
r   = x(1);
beq = x(2);

w	= (1-alffa) * (alffa/(r+delta))^(alffa/(1-alffa));

%% Solve the problem of OLG households given prices, taxes and transfers
% Calculate the sum of lifetime discount factors
surv = cumprod(1-mort);
disc = 1;
J   =  80;

for j = 2:J
    disc = disc + betta^(j-1) * surv(j-1);
end

% Calculate the NPV of lifetime income
lti = 0;
for j = 1:J
    lti = lti + ret*(w + beq)/(1+r)^(j-1);
end

% Find c1 that satisfies the lifetime budget constraint
c1 = lti/disc;

a_0 = 0;

c_all = zeros(J, 1);
h_all = zeros(J, 1);
a_all = zeros(J, 1);

c_all(1) = c1(1);
h_all(1) = 1;
a_all(1) = w + beq + (1+r)*a_0 - c_all(1);

for j = 2:J
    c_all(j) = betta*(1-mort(j-1))*c_all(j-1)*(1+r);
    h_all(j) = ret(j-1);
    a_all(j) = ret(j-1)*w + beq + (1+r)*a_all(j-1) - c_all(j);
end


c = sum(c_all.*N_rel_all)/N_rel;
h = sum(h_all.*N_rel_all)/N_rel;
a = sum(a_all.*N_rel_all)/N_rel/(1+n);

%% Remaining variables
k	= (alffa/(r+delta))^(1/(1-alffa))*h;
inv	= (delta+n)*k;
y	= k^alffa*h^(1-alffa);

F = [a - k;
     beq - sum((N_rel_all(1:end-1)/(1+n) - N_rel_all(2:end)).*a_all(1:end-1))/N_rel*(1+r)];
end

