g=50;n=10;k=.0333;a1=k-1;a2=2-2*k;a3=k;
f=zeros(1,n-1);
v2=[a1 a2 a3];f0=a1*g+a2*g^2+a3*g^3
for l=2:n
v=1:1:2*l-1;
v1=1:1:2*l+1;
va1=v*a1;
va3=v*a3;
ba2=0:-1:2-2*l;
va2=v*a2+ba2;
c1=zeros(2*l-1,2);
c2=zeros(2*l-1,1);
A1=[diag(va1) c1];
A2=diag(va2,1);A21=A2(1:2*l-1,:);A22=[A21 c2]; 
A3=diag(va3,2);A31=A3(1:2*l-1,:);
A4=A1+A22+A31;
e=g .^v1;
v2=v2*A4;
h=v2*e';
gl1=sum(h);
f(l-1)=gl1;
end
v0=[g f0]
f1=[v0 f]
t=1:n;f2=f1(1,1:n);
semilogy(t,f2)
%stem(f2)
for i=1:n
vas(i)=log10(f2(i));
end
stem (vas)
xlabel('k')
ylabel('lg(E(g^k))')

title('The moments of the gain g, G=50, k=0.0333')
