%This script assumes that s, phi, T and r are already defined.
%If s, phi are rows switch them to columns.
if dot([1,-1],size(s))<0 
    s=s';
end
if dot([1,-1],size(phi))<0 
    phi=phi';
end
%
%Calculate the array of s(i+1)-s(i-1) and list of Bank account values.
n=length(s);
dds=s(3:n)-s(1:n-2);
B=(1+r).^(0:T);
%
%Calculate the martingale probabilities.
q=[NaN;((1+r)*s(2:n-1)-s(1:n-2))./dds;NaN];
%
%Initialize v, alpha, beta
v=zeros(n,T+1);
alpha=v; beta=v;
v(:,T+1)=phi;
%
%Fill in v(1,:) and v(n,:) using linear approximation formula:
v(n,1:T)=s(n)*(phi(n)-phi(n-1))/(s(n)-s(n-1))+B(1:T)*(phi(n-1)*s(n)-phi(n)*s(n-1))/(B(T+1)*(s(n)-s(n-1)));
v(1,1:T)=s(1)*(phi(2)-phi(1))/(s(2)-s(1))+B(1:T)*(phi(1)*s(2)-phi(2)*s(1))/(B(T+1)*(s(2)-s(1)));
%
%The iteration.
for k=T:-1:1
    v(2:n-1,k)=(q(2:n-1).*v(3:n,k+1)+(1-q(2:n-1)).*v(1:n-2,k+1))/(1+r);
    alpha(2:n-1,k)=(v(3:n,k+1)-v(1:n-2,k+1))./dds;
    beta(2:n-1,k)=(s(3:n).*v(1:n-2,k+1)-s(1:n-2).*v(3:n,k+1))./(B(k+1)*dds);
    %
    %Portfolio values for s(1) and s(n) as absorbing states
    alpha([1;n],k)=((1+r)*v([1;n],k)-v([1;n],k+1))./(r*s([1;n]));
    beta([1;n],k)=(v([1;n],k+1)-v([1;n],k))/(r*B(k));
end
%
%Display results, flipped top to bottom
disp('s and q:')
disp(flipud([s,q]))
disp('v:')
disp(flipud(v))
disp('alpha:')
disp(flipud(alpha))
disp('beta:')
disp(flipud(beta))