The Jacobian matrix evaluated at the steady state contains elements that are not real or are infinite

I have entered the analytical steady state in the steady state block.

There is only one residual, the last equation. I suspect something is wrong with my bond market clearing equation.

var c cstar n nstar k kstar z zstar r rstar w wstar y ystar b bstar mu mustar q;
predetermined_variables k kstar b bstar;
varexo e estar;
parameters beta theta delta rho rhostar sigma a debtlimit;

beta = 0.99;
delta = 0.025;
theta = 0.36;
a = 1;
rho = 0.906;
rhostar = 0.088;
sigma   = .00852;
debtlimit = 1;

model;
/*
Home Economy
*/
y = exp(z) * k^theta * n^(1 - theta);

c + k(+1) + q * b(+1) = w * n + r * k + (1 - delta) * k + b;
b(+1) = - debtlimit * y;

a/(1-n) = (1/c) * w;
(1/c) = beta * ((1/c(+1)) * (r(+1) + 1 - delta));
beta * (1/c(+1)) = q * (1/c) + mu;

z = rho * z(-1) + rhostar * zstar(-1) + e;
r = theta * exp(z) * k^(theta - 1) * n^(1-theta);
w = (1 - theta) * exp(z) * k^(theta) * n^(-theta);

/*
Foreign Economy
*/
ystar = exp(zstar) * kstar^theta * nstar^(1 - theta);

cstar + kstar(+1) + q * bstar(+1) = wstar * nstar + rstar * kstar + (1 - delta) * kstar + bstar;
bstar(+1) = - debtlimit * ystar;

a/(1-nstar) = (1/cstar) * wstar;
(1/cstar) = beta * ((1/cstar(+1)) * (rstar(+1) + 1 - delta));
beta * (1/cstar(+1)) = q * (1/cstar) + mustar;

zstar = rho * zstar(-1) + rhostar * z(-1) + estar;
rstar = theta * exp(zstar) * kstar^(theta - 1) * nstar^(1-theta);
wstar = (1 - theta) * exp(zstar) * kstar^(theta) * nstar^(-theta);

b(+1) + bstar(+1) = 0;

//y + ystar = c + cstar + k(+1) - (1-delta)*k + kstar(+1) - (1-delta)*kstar;
end;

steady_state_model;
r = (1/beta) + delta - 1;
n = 1 / (1 + (a/(1-theta)) * (1 - delta * theta * (beta/(1 + beta * delta - beta))));
k = (r / (theta * n^(1 - theta)))^(1/(theta-1));
y = k^theta * n^(1 - theta);
w = (1 - theta) * (y / n);
b = -debtlimit * y;
c = y - delta * k;
q = (w * n + r * k - delta * k + b  - c) / b;
mu = (beta/c) - (q/c);


rstar = (1/beta) + delta - 1;
nstar = 1 / (1 + (a/(1-theta)) * (1 - delta * theta * (beta/(1 + beta * delta - beta))));
kstar = (rstar / (theta * nstar^(1 - theta)))^(1/(theta-1));
ystar = kstar^theta * nstar^(1 - theta);
wstar = (1 - theta) * (ystar / nstar);
bstar = -debtlimit * ystar;
cstar = ystar - delta * kstar;
mustar = (beta/cstar) - (q/cstar);

z = 0;
zstar = 0;
end;

model_diagnostics;
resid;
steady;