% DEFINE THE GRID FOR K

kmin = 1; % The minimum value for k
kmax = 20; % The maximum value for k

nk_index = 25:5:200;

for ii = 1:length(nk_index);
    
nk = nk_index(ii); % The total number of grid points
kstep = (kmax - kmin)/(nk-1); %Define the step length of the grid
K = kmin:kstep:kmax; % Define the grid for K

% SET PARAMETERS
alpha = .33; % A production function parameter
beta = .95; % The Discount factor
delta = 0; % The rate of depreciation
tol = 1e-5; % The convergence tolerance
maxit = 1000; % The maximum number of iterations


% Initialize miss difference and iteration counter
dd = 1;
it = 1;

% Provide initial guess for gi, the policy function index
gi = 1:nk;
gi = gi';

% Define R for single value function iteration

KS = repmat(K',[1 nk]);
KC = repmat(K',[1 nk])';

%%% NOTE: THIS IS WHERE DELTA MUST BE ADDED!!!
R = KS.^alpha + (1-delta)*KS - KC;


while (it <= maxit)
    J = zeros(nk,nk);
    r = zeros(nk,1);
for i=1:nk;
   J(i,gi(i)) = 1;
   r(i) = R(i,gi(i));
end;


V = (eye(nk) - beta*J)\r;
OBJ = R + beta*ones(nk,1)*V';
[Vnew,ginew] = max(OBJ,[],2); % maximize the objective function
dd = max(abs(gi-ginew)); % Calculate the miss distance
gi = ginew;

if (max(gi) == nk)
    disp('Warning: Potential Corner Solution.  Increase kmax');
end;

if (min(gi) == 1 && it>1)
    disp('Warning: Potential Corner Solution.  Decrease kmin');
end;
if (dd<tol);
    break;
end


it = it+1;

end;

KSTAR(ii) = K(gi(1));

end