function theta = dgrowc(kmin,kmax,n,param);
%DGROWC  This program computes the collocation approximation for the
%        consumption function of the deterministic growth model.
%        theta = dgrowc(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);

% Call a subroutine for quadrature abscissas and weights.

x     = linspace(kmin,kmax,n)';
m     = n;
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   = (beta*alpha*lambda)^(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;

% Use Newton's method to solve the n residual equations for theta.

for j=1:100;
  C     = K*theta;
  xp    = lambda*x.^alpha-C;
  Kp    = (xp*ones(1,n)).^(ones(m,1)*[1:n]);
  Cp    = Kp*theta;
  F     = Cp-beta*alpha*lambda*xp.^(alpha-1).*C;
  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;
    Kp    = (xp*ones(1,n)).^(ones(m,1)*[1:n]);
    Cp    = Kp*theta;
    Fp    = Cp-beta*alpha*lambda*xp.^(alpha-1).*C;
    dF(:,i) = (Fp-F)/del(i,i);
  end;
  theta = thold - dF\F;
  if norm(theta-thold)<1e-10; break; end;
end;
