[X]=e633

function [X]=e633
%
%  Example 6.3.3
%
pt=[.05 .05 .10 .15 .50 .10 .05];
pd=[.2 .5 .2 .1]; h=1; K=25; c=4; pen=20; N=6; NS=N+1;
n=length(pt); t=cumsum(ones(1,n)); ET=t*pt';
m=length(pd); Li=[];
for i=0:m-1
   suma=0;
   i1=max(i,1);
   for k=m:-1:i1
      suma=suma+(k-i)*pd(k);
   end
   suma=suma*pen; Li=[Li suma];
end
cost=[]; sx=0:N; X=zeros(NS,2); X(:,1)=sx';
for s=0:N
P=zeros(NS,NS); tau=[]; r=[];
for i=s:N
   for j=1:i
      ij=i-j;
      if ij > 0 
      if ij < 5
      P(i+1,j+1)=pd(ij);
      end
      end
   end
tau=[tau ET]; r1=i*h*ET;
if i < m
r1=r1+Li(i+1);
end
P(i+1,1)=1-sum(P(i+1,:)); r=[r r1];
end
for i=s-1:-1:0
   P(i+1,NS)=1; tau=[0 tau]; r1=K+c*(N-i); r=[r1 r];
end
IP=eye(NS)-P; IP(:,1)=ones(NS,1); IP=inv(IP); pi=IP(1,:);
up=pi*r'; down=pi*tau'; 
g=up/down;
cost=[cost g];
end
X(:,2)=cost';