










%hold off
t0=clock;
k1=.5;g=50;n3=3;
a=(1+k1*(g-1))/(g);
n1=600;%p1=2.756250000000000e-001
 %5.677930563756942e-002;% (g=80)  % 0.05145197461537;
m=(1-k1)*(1-a);
n=a*(1-k1);
b=k1*(1-a);
p1=a^(1/(1-k1))
p2=p1^2*(m+b)/n;
v=[p1 p2];x=zeros(1,n1-2);
p=[v x];w=zeros(1,n1);w(n3)=a^(n3/(1-k1));q=n3/((n3+1)*n);
v1=a^(n3/(1-k1));s=a^(k1/(1-k1));t=(1-a)/(1-k1);v2=n3*t*v1*s;v3=.5*v2*t*s*(n3+1+k1);
v4=[v1 v2 v3];
for i=3:n1,
suma=0;
k=fix(i/2);
l=1;
while (l<k)
suma=suma+p(l)*p(i-l);
l=l+1;
end
if i-2*k>0
	suma=suma+p(k)*p(i-k);
	p(i)=suma*(2*m+i*b)/((i-1)*n);
else
	suma=suma+p(k)^2/2;
	p(i)=suma*(2*m+i*b)/((i-1)*n);
end
end
duratad=etime(clock,t0)
resd=p;
%stem(resd)
x=(1:n1);
whitebg;
%plot(x,resd)
loglog(x,resd,'k')
xlabel('m');
ylabel('P(m/1,P(m/3)');
hold on
pack
%sigma=0;
for l=n3+1:n1
sigma=0;%w(n3)=a^(n3/(1-k1));
	for j=1:l-n3
	sigma=sigma+p(j)*w(l-j);
	end
  w(l)=sigma*q*(l*b+m*(n3+1))/(l-n3);
end

loglog(x,w,'k');
grid	
title('The probabilities P(m/1) and P(m/3),k=0.5, G=50  ')
%white bg;
cc=conv(p,p);
ee=conv(cc,cc);
ff=conv(ee,ee);
dd=w(1,8:200)-ff(1,1:193);
duratae=etime(clock,t0)
mc=w(1,n3:n3+15);
%inverthardcopy on