MATLABで"添字インデックスは、実数の正の整数か、論理値のいずれかでなければなりません。"というエラーの対処方法をご教示頂けないでしょうか。

###前提・実現したいこと
マクロ経済学のモデルを構築し、そのモデルにおける定常状態を求めたいと考えております。

###発生している問題・エラーメッセージ
"添字インデックスは、実数の正の整数か、論理値のいずれかでなければなりません。"

###該当のソースコード

%%%%パラメーター%%%%%%%%%%%
alpha=0.4; %%資本分配率
sigma=1/1.5;
beta=0.99;
delta=0.05; %%資本減耗率 *3期間に合わせる
A=1;%%経済の規模
e=10^(-8); %%収束条件
eta=0.8; %%定常状態の計算の際に新たな値を見つけるときの
KI=1; %% kの初期値
t=100; %%期間
gamma=0.5; %% gamma=1のとき対数効用と一致
nn=0.02;%%人口成長率
ttyp=0.09;%%若年期保険料支払い率 17.828%の半分の8.914%
ttmp=0.09;%%中年期保険料支払い率
wtax=0.2; %%賃金税
ctax=0.08; %%消費税
gp=0.15;%%政府はGDPの一定割合15%を支出
m=100; %Dのpathの試行回数
sp=0.5;%政府が年金を補う比率
bb=0.5;%年金の所得代替率
ee=2;%若年勤労期の労働供給量を1としたときの中年勤労期の労働供給量
ph=1.2;%middleと比べてph倍の労働を供給する
pm=1;%%基準であるmiddleの供給量
pl=0.8;%%%middleと比べてpl倍の労働を供給する

%%%%%%%%%% find the steady state%%%%%

ks=KI; %%kの初期値
for i=1:200 %% if(ksと計算されたkaの差がeより大きかったら) 新たなksを使って200回まで計算し直す、 ksと計算されたkaがeより小さかったらbreak

%%企業部門
ys=A*ks^(alpha);
rs=alpha*A*ks^(alpha-1)-delta;
ws=(1-alpha)*A*ks^alpha;

%%年金部門
bsH=((((1+nn)^2)*ttyp*ws*ph+(1+nn)*ttmp*ws*ee*ph)*bb)/(2*(1-sp)); %%Highの年金受け取り
bsM=((((1+nn)^2)*ttyp*ws*pm+(1+nn)*ttmp*ws*ee*pm)*bb)/(2*(1-sp)); %%Middleの年金受け取り
bsL=((((1+nn)^2)*ttyp*ws*pl+(1+nn)*ttmp*ws*ee*pl)*bb)/(2*(1-sp)); %%Lowの年金受け取り

%%High skill家計部門
%young
cityH=((1-wtax-ttyp)*ws*ph+((1-wtax-ttmp)*ee*ws*ph)/(1+rs)+(bsH/((1+rs)^2)))/(1+ctax+(beta^(1/gamma))*((1+rs)^((1-gamma)/gamma))*(1+ctax)+(beta^(2/gamma))*((1+rs)^((2-2*gamma)/gamma))*(1+ctax));
sityH=(1-wtax-ttyp)*ws*ph-(1+ctax)*cityH;
%middle
ci1tmH=((1-wtax-ttmp)*ee*ws*ph+(1+rs)*sityH+bsH/(1+rs))/(1+ctax+(beta^(1/gamma))*((1+rs)^((1-gamma)/gamma))*(1+ctax));
si1tmH=(1-wtax-ttmp)*ee*ws*ph+(1+rs)*sityH-(1+ctax)*ci1tmH;
%%old
ci2toH=(1+rs)*si1tmH/(1+ctax)+bsH/(1+ctax);

%%Middle skill家計部門
%young
cityM=((1-wtax-ttyp)*ws*pm+((1-wtax-ttmp)*ee*ws*pm)/(1+rs)+(bsM/((1+rs)^2)))/(1+ctax+(beta^(1/gamma))*((1+rs)^((1-gamma)/gamma))*(1+ctax)+(beta^(2/gamma))*((1+rs)^((2-2*gamma)/gamma))*(1+ctax));
sityM=(1-wtax-ttyp)*ws*pm-(1+ctax)*cityM;
%middle
ci1tmM=((1-wtax-ttmp)*ee*ws*pm+(1+rs)*sityM+bsM/(1+rs))/(1+ctax+(beta^(1/gamma))*((1+rs)^((1-gamma)/gamma))*(1+ctax));
si1tmM=(1-wtax-ttmp)*ee*ws*pm+(1+rs)*sityM-(1+ctax)*ci1tmM;
%%old
ci2toM=(1+rs)*si1tmM/(1+ctax)+bsM/(1+ctax);

%%Low skill家計部門
%young
cityL=((1-wtax-ttyp)*ws*pl+((1-wtax-ttmp)*ee*ws*pl)/(1+rs)+(bsL/((1+rs)^2)))/(1+ctax+(beta^(1/gamma))*((1+rs)^((1-gamma)/gamma))*(1+ctax)+(beta^(2/gamma))*((1+rs)^((2-2*gamma)/gamma))*(1+ctax));
sityL=(1-wtax-ttyp)*ws*pl-(1+ctax)*cityL;
%middle
ci1tmL=((1-wtax-ttmp)*ee*ws*pl+(1+rs)*sityL+bsL/(1+rs))/(1+ctax+(beta^(1/gamma))*((1+rs)^((1-gamma)/gamma))*(1+ctax));
si1tmL=(1-wtax-ttmp)*ee*ws*pl+(1+rs)*sityL-(1+ctax)*ci1tmL;
%%old
ci2toL=(1+rs)*si1tmL/(1+ctax)+bsL/(1+ctax);

%%政府部門
ds=-(gp*ys+((sp*bb*ttyp*ws)*(ph+(1+nn)*ee*ph+pm+(1+nn)*ee*pm+pl+(1+nn)*ee*pl))/(2*3*(2+nn))-wtax*ws*(((1+nn)*(ph+pm+pl)+ee*(ph+pm+pl))/(3*(2+nn)))-ctax(((1+nn)/(3*(2+nn)))*(cityH+cityM+cityL)+(1/(3*(2+nn)))*(ci1tmH+ci1tmM+ci1tmL)+(1/(3*(1+nn)*(2+nn)))*(ci2toH+ci2toM+ci2toL)))/rs;

%%市場均衡
ka=(1/(3*(2+nn)))*(sityH+sityM+sityL)+(1/(3*(2+nn)*(1+nn)))*(si1tmH+si1tmM+si1tmL)-ds;

if abs(ka-ks)<e
    break
end
ks=eta*ks+(1-eta)*ka;