function [X,u,mu,A,B]=nlsteady(beta,D,E,param,param3,crit,damp,maxit,perdev, ..
                           X0,u0,mu0);
%NLSTEADY computes the steady state values X, u, mu of the following problem:
%            
%       (P)      max sum{beta^t u(x(t),param) }   param=[param1;param2]
%                s.t.  x(t) = D*X(t)+E*u(t)       X(t): ns x 1,  u(t): nc x 1
%                      X(t+1) = h(X(t),u(t),param3)    
%                      X(1,t)=1  for all t
%                      mu(t)=lagrange mult. for X(t+1) = h(X(t),u(t))
%                            (mu(1)=-9999)
%                                                      
%  INPUTS:
%     1. parameters for problem (P): beta,D,E,param,param3
%     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. perdev=the positive vector of length ns+nc with elements equal to the 
%        percent deviations desired in approximating numerical derivatives.
%     6. 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.  
%     
%        Caution: If there are any constraints (with the exception of the first
%                 one which are of the form: X(i,t+1)=X(i,t), the algorithm used
%                 will break down.  Reformulate the problem setting X(i) 
%                 constant wherever it appears.
%
%  OUTPUTS:
%     1. X,u,mu: steady state values of the problem (P)
%
[row,ns]=size(D);
[row,nc]=size(E);
%
%  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;
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) = h(X(t),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*(dh(X(t+1),u(t+1))/dXs(t+1))*mu(t+1) =0
%  u:    E'*uj -(dh(X(t+1),u(t+1))/du(t+1))*mu(t) =0
%
for i=1:maxit;
  adm=diag(.01*perdev.*[1;X0;u0]);
  arg=find(abs([1;X0;u0])<.0001);
  for j=1:length(arg);
    adm(arg(j),arg(j))=.001;
  end;
  h=hfn([1;X0;u0],param3);
  for j=1:ns+nc;
    jac(j,:)=(hfn([1;X0;u0]+adm(:,j),param3)'-hfn([1;X0;u0]-adm(:,j), ..
              param3)')/(2*adm(j,j));
    H(j,j)=((hfn([1;X0;u0]+adm(:,j),param3)'-2*h'+hfn([1;X0;u0]-  ..
              adm(:,j),param3)')/(adm(j,j)^2))*[1;mu0];
    for k=1:j-1;
      H(j,k)=((hfn([1;X0;u0]+adm(:,j)+adm(:,k),param3)'-hfn([1;X0;u0]+ ..
       adm(:,j)-adm(:,k),param3)'-hfn([1;X0;u0]-adm(:,j)+adm(:,k),param3)'+ ..
       hfn([1;X0;u0]-adm(:,j)-adm(:,k),param3)')/(4*adm(j,j)*adm(k,k)))*[1;mu0];
      H(k,j)=H(j,k);
    end;
  end;
  [us,uj,uh]=ufn(D*[1;X0]+E*u0,param); 
  M1=[eye(nss)-jac(2:ns,2:ns)'          -jac(ns+1:ns+nc,2:ns)'      zeros(nss);
      beta*C1'*D'*uh*D*C1-beta*H(2:ns,2:ns),    ...
          beta*C1'*D'*uh*E-beta*H(2:ns,ns+1:ns+nc) eye(nss)-beta*jac(2:ns,2:ns);
      E'*uh*D*C1-H(ns+1:ns+nc,2:ns),   ...
          E'*uh*E-H(ns+1:ns+nc,ns+1:ns+nc)     -jac(ns+1:ns+nc,2:ns)];
  M2=[-X0+h(2:ns);
      -(eye(nss)-beta*jac(2:ns,2:ns))*mu0-beta*C1'*D'*uj;
      jac(ns+1:ns+nc,2:ns)*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 NLSTEADY.'); 
else;
  sprintf('Convergence achieved in NLSTEADY at iteration: %g',i)
end;
X=[1;X1];
u=u1;
mu=[-9999;mu1];
%
%  The first order Taylor expansion of h(X(t),u(t)) results in A*X(t)+B*u(t),
%  where A and B are given by:
%
B=jac(ns+1:ns+nc,:)';
A0=X-jac(1:ns,:)'*X-B*u;
A=[         1            zeros(1,nss);
   A0(2:ns)+jac(1,2:ns)' jac(2:ns,2:ns)'];
