a=1.05;b=0.01;q=1;r=1;x0=10;sN=5;N=500;

sinfty = roots([b^2 (r-a^2*r-b^2*q) -q*r]), 
sinfty = sinfty(1);

kinfty = a*b*sinfty/(r+b^2*sinfty)
aclosed = a-b*kinfty

[x,u,K,S]=scaopt(a,b,q,r,sN,x0,N);

k=1:N;
subplot(1,3,1)
plot(k,S(1:end-1),k,repmat(sinfty,1,length(k)));
legend('Optimal','Steady-state optimal')
xlabel('k')
ylabel('s_k')

subplot(1,3,2)
plot(k,K,k,repmat(kinfty,1,length(k)))
legend('Optimal','Steady-state optimal')
xlabel('k')
ylabel('K_k')

G = ss(a,b,1,0,'Ts',-1)
Gclosed = feedback(G,kinfty)

subplot(1,3,3)
[yinfty,k,xinfty]=initial(Gclosed,x0,k);
plot(k,x(1:end-1),k,xinfty)
legend('Optimal','Steady-state optimal')
xlabel('k')
ylabel('x_k')