gridpoints=700;
gridend=3.5;
step_size=0.5;
mu=0.3;
cv=0.3;
r=0.06;
eps=2e-6;
vardelta=0.1;
varbeta=0.4;
varpsi=(1-vardelta)/(1+r);
varsigma=.1;
xi=(r+vardelta)/varbeta;
grid=-gridend:(2*gridend/(gridpoints-1)):gridend;
diff=ones(length(grid),1);
profits=max(0,xi*(exp(varbeta*grid)-varbeta*exp(grid)));
valfold=profits;
density=normpdf(grid, log(1-vardelta) , varsigma );
density=density/sum(density);
expxy=exp(-grid);
half=ceil(length(grid)/2);
while max(abs(diff))>eps;
  omega=1/xi*exp(-varbeta*grid).*( max(valfold) - valfold );
  summation=zeros(1,length(grid)); 
  for y=1:length(grid);
     ints(y)=sum((gamcdf(0.00:0.001:omega(y), cv^(-2), mu*cv^2))/max(omega(y),.001))/1000;
     c1=circshift(fliplr(density)',[half+y-1,0])';
     c2=circshift(fliplr(expxy'),[y+half,0])';
     if y>half 
        c1(1:(y-half))=0;
     end
     if y<half 
        c1((half+y):end)=0;
     end
     summation=summation+c1.*c2.*(valfold(y) + xi*exp(varbeta*grid(y))*ints(y));
  end
  % Update the value function only slowly; Because the value function is convex in some parts, the value function iteration can get stuck in cycles if step size =1 (as you would normally do) 
  valf=(1-step_size)*valfold+step_size*(profits+varpsi*summation);
  diff=valf-valfold;
  valfold=valf;
end;
