function [X,u,mu]=steady(beta,A,B,D,E,param,crit,damp,maxit,X0,u0,mu0);
%STEADY computes the steady state values X, u, mu of the following problem:
%            
%       (P)      max sum{beta^t u(x(t),param) }
%                s.t.  x(t)=D*X(t)+E*u(t)       X(t): ns x 1,  u(t): nc x 1
%                      X(t+1)=A*X(t)+B*u(t)    
%                      X(1,t)=1  for all t
%                      mu(t)=lagrange mult. for X(t+1)=A*X(t)+B*u(t)
%                      NOTE: only mu elements for non-redundent or nontrivial
%                            constraints can be computed.  The others will 
%                            be set at -9999.
%                                                      
%  INPUTS:
%     1. parameters for problem (P): beta,A,B,D,E,param
%     2. crit=convergence criterion
%     3. damp=damping factor for slower (but more stable) convergence. damp=1
%        is full damping (no changes allowed), damp=0 is no damping.
%     4. maxit=maximum # of iterations
%     5. X0,u0,mu0: initial conditions for steady states.  These are optional
%        arguments.  If you do not have a guess for any subset of these 
%        vectors, put a -9999 in the argument.  If length(mu0)=number of 
%        states then the elements of mu0 corresponding to trivial or redundent
%        constraints are ignored. If length(mu0)=number of independent
%        constraints, all are used in computation.
%
%  OUTPUTS:
%     1. X,u,mu: steady state values of the problem (P)
%
[ns,col]=size(A);
[row,nc]=size(B);
%
%  Since I assume the system includes a constant term (i.e. X(1,t)=1 for all 
%  t) and thus a constraint:  X(1,t+1)=X(1,t), I reduce the system to avoid
%  inverting a singular matrix:
%
if X0==-9999;  X0=ones(ns-1,1); else; X0=X0(:); X0=X0(2:ns,1); end;
if u0==-9999;  u0=ones(nc,1); else; u0=u0(:); end;
if mu0==-9999; mu0=ones(ns-1,1); else; mu0=mu0(:); mu0=mu0(2:ns,1); end;
A0=A(2:ns,1);
A=reducmat(A,1,1);
B=reducmat(B,1,0);
C1=reducmat(eye(ns),0,1);
nss=ns-1;
tnn=2*nss+nc;
e=ones(tnn,1);
M20=9999999*ones(tnn,1);
%
%  Newton-Raphson is used to solve the first-order conditions of 
%     
%      max( sum{beta^t*ufn(D*(C0+C1*Xs(t))+E*u(t)) } )
%      st.  Xs(t+1)=A0 + As*Xs(t)+Bs*u(t)  (mult: mu)   
%           X(t)= C0+C1*Xs(t)
%
%  where subscript "s" refers to the system without the first constraint;
%  Note: in the program, the reduced matrices are not renamed to save memory.
%
%  These conditions (in addition to the constraints) are:
%      
%  Xs:   beta*C1'*D'*uj(t) +mu(t) -beta*As'*mu(t+1) =0
%  u:    E'*uj -Bs'*mu(t) =0
%
if rank(eye(nss)-A)>=nss; 
  X0=(eye(nss)-A)\(A0+B*u0);
end;
[us,uj,uh]=ufn(D*[1;X0]+E*u0,param);
mu0=-(eye(nss)-beta*A')\(beta*C1'*D'*uj);
for i=1:maxit;
  [us,uj,uh]=ufn(D*[1;X0]+E*u0,param); 
  M1=[eye(nss)-A               -B              zeros(nss);
   beta*C1'*D'*uh*D*C1   beta*C1'*D'*uh*E    eye(nss)-beta*A';
   -E'*uh*D*C1            -E'*uh*E              B'];
  M2=[A0-(eye(nss)-A)*X0+B*u0;
      -(eye(nss)-beta*A')*mu0-beta*C1'*D'*uj;
       -B'*mu0+E'*uj];
  %  Are we diverging?
  if (abs(M2) > abs(M20)); 
    disp('The procedure is diverging.  Try new starting values.')
    return
  end;
  MIM=(1-damp)*(M1\M2);
  X1=X0+MIM(1:nss,1);
  u1=u0+MIM(ns:nss+nc,1);
  mu1=mu0+MIM(ns+nc:tnn,1);
  %  check for convergence:
  xum0=[X0;u0;mu0]
  if max(abs( ([X1;u1;mu1]-xum0)./max(abs(xum0),e) )) < crit; break; end;
  X0=X1;
  u0=u1;
  mu0=mu1;
  M20=M2;
end;
if i>=maxit; 
  disp('Iteration limit reached in STEADY.'); 
else;
  sprintf('Convergence achieved in STEADY at iteration: %g',i)
end;
X=[1;X1];
u=u1;
mu=[-9999;mu1];
