function theta = dgrowl(kmin,kmax,n,param);
%DGROWL  This program computes the Least Squares approximation for the
%        consumption function of the deterministic growth model.
%        theta = dgrowl(kmin,kmax,n,param) returns the vector theta with
%        the minimum value of $k$ (kmin), the maximum value of $k$ (kmax), 
%        the degree of the polynomial (n), and the parameters of the 
%        problem (param) as inputs.
%

%        Ellen McGrattan, 8-30-96

alpha = param(1);
beta  = param(2);
lambda= param(3);
m     = param(4);
b     = beta*alpha*lambda;

% Call a subroutine for quadrature abscissas and weights.

[x,w] = qgausl(kmin,kmax,m);
K     = (x*ones(1,n)).^(ones(m,1)*[1:n]);

% Make an initial guess for c(k) and translate it to a guess for theta:

kss   = b^(1/(1-alpha));
if (kss<kmin | kss>kmax);
  disp(['WARNING: Steady state capital not in [', ...
        num2str(kmin),',',num2str(kmax),']'])
end;
css   = lambda*kss^alpha-kss;
c     = css * sqrt(x/kss);
c     = (1-alpha*beta)*lambda*x.^alpha;
theta = K\c;

% Evaluate the n residual equations.

for j=1:100;
  C     = K*theta;
  xp    = lambda*x.^alpha-C;
  Kp    = (xp*ones(1,n)).^(ones(m,1)*[1:n]);
  dKp   = ([1:m]'*ones(1,n)) .* (xp*ones(1,n)).^(ones(m,1)*[0:n-1]);
  Cp    = Kp*theta;
  R     = Cp-b*xp.^(alpha-1).*C;
  tem   = dKp*theta - b*((alpha-1)*xp.^(alpha-2).*C+xp.^(alpha-1));
  dR    = K.*(tem*ones(1,n))+Kp;
  F     = sum( dR .* ( (w.*R) * ones(1,n)) )';
  thold = theta;
  del   = diag(max(1e-8,abs(theta)*.0001));
  for i=1:n;
    theta    = thold;
    theta(i) = theta(i)+del(i,i);
    C     = K*theta;
    xp    = lambda*x.^alpha-C;
    Cp    = (xp*ones(1,n)).^(ones(m,1)*[1:n])*theta;
    R     = Cp-b*xp.^(alpha-1).*C;
    tem   = dKp*theta - b*((alpha-1)*xp.^(alpha-2).*C+xp.^(alpha-1));
    dR    = K.*(tem*ones(1,n))+Kp;
    Fp    = sum( dR .* ( (w.*R) * ones(1,n)) )';
    dF(:,i) = (Fp-F)/del(i,i);
  end;
  theta = thold - dF\F;
  if norm(theta-thold)<1e-10; break; end;
end;


