"Array indices must be positive integers or logical values." while using golden section search with an interpolation?

The function starts going wrong as I hit line 50, the one that starts with w(i,j). What am i doing wrong? A seperate function is implemented in the code called gss which represents the golden search function, code posted below:

%This script performs value function iteration for the baseline model.

clear all
global delta eta alpha beta
delta=1; %depreciation rate
eta=1; %utility parameter
alpha=.27; %production function parameter
beta=.994; %patience parameter
eps=.01; %partial tolerance parameter
n=250; %grid size
kstar=((1-beta*(1-delta))/(alpha*beta))^(1/(alpha-1)); %steady-state capital
Kl=.75*kstar; %lower bound for capital grid
Ku=1.25*kstar; %upper bound for capital grid
kstep=(Ku-Kl)/(n-1); %step size for capital grid
kgrid=Kl:kstep:Ku; %grid for capital stock
v0=zeros(n,1); %initialize the value function
dist=1; %initialize the distance for the value function
it=1; %initialize the counter
tic %initialize the clock
eps2=.000001;


while (dist>eps*(1-beta))
    for i=1:length(kgrid)
        %Step 1
        for j=1:length(kgrid)
            ct=kgrid(i)^alpha+(1-delta)*kgrid(i)-kgrid(j); %consumption
            %if i=1, kgrid(i)=0.1237. If I choose a kgrid of j=1, the next
            %period will also be 0.1237. As we mvoe to the next period,
            %kgrid(j)=0.1246. 
            
            ct=ct*(ct>0)+eps2; %correct for potential negative consumption
            
            %ij=[kgrid(j-1) krid(j+1)]- needed component; 
            %yi_given=interpl(x,y,xi_given,'linear'); x&y are given strings
            %of values or you can represent it as below 
            %yp=interpl(x,y,xp,'linear'); xp=[0:0.25:10]; The interpolation
            %is a way to make the points finer to match the now finer gird
            %neeeded for golden section searhc as the gss method demand
            %the input of an interval with only upper and lower bound with
            %no step size.

            %you set i=1 and go through all these js. For every single i,
            %you will try to find the j that maxmizes that value function.
            %Essentially, you keep the i fixed, i=1 and go through all the
            %columns represneted by j. For instance, it goes from (1,1) to
            % (1,2 & n) for i=1. You are trying to find the max within that 
            %row through step 2.1. 
            
            w(i,j)=gss(utilchp4(ct,eta)+beta*interpl(1:j,v0(j,1),kgrid(j-1):krid(j+1),'linear'),[kgrid(j-1) krid(j+1)]); 
            
            %value function are 
            %monotonically increasing. This function stores the value 
            %created by capital as indicated by i and j.utilchp4(ct,eta)
            %give you a instaneous return. Use golden section search to
            %lienar interpolate that entire right hand side. Professor
            %highlighted v0(j,1) when speaking about linear interpolation.
            
        end
        %Step 2
        
        %find the max of w for every i. When you find that, that points to
        %jstar. The location of the column of the biggest W essentially
        %represent the location of the column of the biggest j. 
        
        jind=find(max(w(i,:))==w(i,:)); %find the j that max rhs of Bellman.
        %when i=2, jstar will be 1, which indicate the first value in the
        %whole string of values as laid out by kgrid. 
        
        vnew(i,1)=w(i,jind); %update the value function for each i. This 
        %tells you where in W that max is. 
        
        %Step 3 
        h(i,1)=jind; % form your policy function. For every i, store jstar.
        %For each iteration, store jstar. Check line 85. 
        
    end
    
    dist=max(abs(v0-vnew)); %calculate distance between value functions
    diststore(it)=dist; %store the distance between value function iterations
    %pause(.2)
    %plot(1:1:it,diststore)
    it=it+1; %update counter
    v0=vnew;%update the value function and repeat. The method discussed in 
    %in class notes as we keep on using the last kprime to get converged
    %toward the kprime 
   
end
kprime=kgrid(h); %given policy pointer h, find the capital stock value
disp('time to completion')
toc 

The golden section search code is provided below (validated with no problem):

function [x]=gss(fun,X)

eps= 0.000001;

p=(sqrt(5)-1)/2; %the golden ratio 
q=1-p;

% Step 1 
   a=X(1); %original lower bound 
   d=X(length(X)); %original upper bound 
   b=p*a+q*d; %lower bound combined with golden ratio calculation 
   c=q*a+p*d;%upper bound combined with golden ratio calculation
   fb=fun(b); 
   fc=fun(c);

% step2+3
   while abs(d-a)>eps*max(1,(abs(a)+abs(c))) %make the interval smaller and smallers 
       
      if fb<fc % choose [b,d] as next interval
         a=b; %cut the range from the right hand side 
         b=c; 
         fb=fc;
         c=p*b+q*d;
         fc=fun(c);
      else % choose [a,c] as next interval  
         d=c; %cut the range from the left hand side 
         c=b;
         fc=fb;
         b=p*c + q*a;
         fb=fun(b);
      end
   end

   %as the function keeps on iterating.... 
   if fb>fc
       x=b; 
   else 
       x=c; 
   end
   
   %disp('To prove that it actually equals to 0')
   %the_function_result=fun(x)
   
end

Answer 1

[Walter Roberson]

            w(i,j)=gss(utilchp4(ct,eta)+beta*interpl(1:j,v0(j,1),kgrid(j-1):krid(j+1),'linear'),[kgrid(j-1) krid(j+1)]);

What are you expecting the kgrid portion of that code to do when j is 1 or j is length(kgrid)?

[Kevin Zhou]

That is a good question. Idealistically, I just want them to start and stop at j is 1 and lenght(kgrid). Would you know how to modify it?