alpha=.35;
beta=.98;
delta=.025;
sigma=2;

gamma1=4.95;
gamma2=5.05;
gamma_bar=5;

g=8; % Dimension of X
eps=1e-4; % Criterion for value function iteration

% ii.
kstar=(alpha*gamma_bar/((1/beta)-1+delta))^(1/(1-alpha))

% iii.
k_min=0.5*kstar;
k_max=1.5*kstar;
X=linspace(k_min,k_max,g);

% iv.
C1=(gamma1*X'.^(alpha)+(1-delta)*X')*ones(1,g)-ones(g,1)*X;
C2=(gamma2*X'.^(alpha)+(1-delta)*X')*ones(1,g)-ones(g,1)*X;

U1=(C1.^(1-sigma)-1)/(1-sigma);
U1(C1<0)=-Inf;

U2=(C2.^(1-sigma)-1)/(1-sigma);
U2(C2<0)=-Inf;

% v.
v1_old=zeros(g,1);
v2_old=zeros(g,1);

v1=max(U1,[],2);
v2=max(U1,[],2);
 
% vi.
i=0;
while max(max(abs((v1-v1_old)./v1)), max(abs((v2-v2_old)./v2)))>eps
    v1_old=v1; v2_old=v2;
    v1=max(U1+beta*.5*ones(g,1)*(v1+v2)',[],2);
    v2=max(U2+beta*.5*ones(g,1)*(v1+v2)',[],2);
    i=i+1;
end
i
% plot(X,[v1 v2])

% vii
[v1,d1]=max(U1+beta*.5*ones(g,1)*(v1+v2)',[],2);
[v2,d2]=max(U2+beta*.5*ones(g,1)*(v1+v2)',[],2);
    
subplot(2,1,1); plot(X,[X(d1); X(d2)]);
legend('\gamma_1','\gamma_2');
xlabel('k'); ylabel('k`');
subplot(2,1,2); plot(X,[X(d1)-X; X(d2)-X]);
legend('\gamma_1','\gamma_2');
xlabel('k'); ylabel('Increase in k');

