function [cmax,v1] = function(param,xa,ya,maxit)
%QR2DP  This file accompanies the article ``Solving the Stochastic Growth
%       Model with a Finite Element Method'' and is used to compute the
%       solution to the dynamic program:
%
%       v(k,ome) = max {c^(1-tau)/(1-tau)-beta E v(k',ome')}
%                   c 
%
%       subject to  c+i  = exp(ome)*k^alpha
%                   k'   = (1-delta)*k+i
%                   ome' = rho*ome + eps    eps ~ N(0,sigma)
%
%       [c,v] = qr2dp(param,xa,ya) returns the optimal consumption levels,
%       c, and the corresponding value function values, v.
%       The vector param is 6 x 1 are the values for alpha, beta, delta,
%       rho, sigma, and tau (in that order).  The vector xa is the partition
%       for the capital stock (k) and the vector ya is the partition for the
%       technology shocks (ome).  Thus, the (i,j) elements of c and v
%       correspond to the ith element of xa and the jth element of ya.
%       [c,v] = qr2dp(param,xa,ya,maxit) allows the user to specify the
%       number of iterations that should be taken.  The default is 1000.

%       Ellen McGrattan, 6-14-93


if nargin<4; maxit=1000; end
alpha   = param(1);
beta    = param(2);
delta   = param(3);
rho     = param(4);
sige    = param(5);
tau     = param(6);
pi      = 3.1415926540;
v0      = (xa.^alpha*ya').^(1-tau)/(1-tau)/(1-beta);
nep     = 10;
[ep,we] = qgausl(nep,-2.88*sige,2.88*sige);
wgt     = exp(-ep.*ep/(2*sige*sige)).*we/sqrt(2*pi*sige*sige);
sum1    = sum(wgt);
wgt     = wgt/sum1;
nx1     = length(xa);
ny1     = length(ya);
nx      = nx1-1;
ny      = ny1-1;

for it = 1:maxit
  for j=1:ny1
    for i=1:nx1
       f   = xa(i)^alpha*ya(j);
       ca  = [.59*f:f/200:f,f];
       for k=1:length(ca)
         % get next period state and integral 
         xt   = (1-delta)*xa(i)+ya(j)*xa(i)^alpha-ca(k);
         i1   = fix( max(xt-xa(1:nx),0)./(xt-xa(1:nx))+.2 );
         i1   = max(1,sum(i1));
         i2   = i1+1;
         dx   = xa(i2)-xa(i1);
         yt   = ya(j)^rho*exp(ep); 
         sum1 = 0;
         for l = 1:nep
           j1   = fix( max(yt(l)-ya(1:ny),0)./(yt(l)-ya(1:ny))+.2 );
           j1   = max(1,sum(j1));
           j2   = j1+1;
           dy   = ya(j2)-ya(j1);
           sum1 = sum1+(v0(i1,j1)*(xa(i2)-xt)*(ya(j2)-yt(l))/dx/dy ...
                       +v0(i2,j1)*(xt-xa(i1))*(ya(j2)-yt(l))/dx/dy ...
                       +v0(i2,j2)*(xt-xa(i1))*(yt(l)-ya(j1))/dx/dy ...
                       +v0(i1,j2)*(xa(i2)-xt)*(yt(l)-ya(j1))/dx/dy)*wgt(l);
         end
         r(k) = ca(k).^(1-tau)/(1-tau)+beta*sum1;
       end
       r         = r(1:length(ca));
       imax      = find(r==max(r));
       imax      = imax(1);
       cmax(i,j) = ca(imax);
       v1(i,j)   = r(imax);
    end
  end
  d = norm(v1-v0);
  if (d<1e-3); break; end
  v0 = v1;
end
disp(sprintf('The norm(v1-v0) at iteration %g was %g.',it,d)) 


