function y=GFP(k1,k2,N,M1,M2);
k=2*k1/k2;
a(1)=1;
s=a(1);
for i=1:N-1
	a(i+1)=a(i)*k*(((i+1)*i)^(-1/2))*((i)^(1/3)+1)^2;
	s=s+a(i+1);
end
for i=1:N
	b(i)=a(i)/s;
end
s1=0;
s2=0;
for i=1:N
	s1=s1+i*b(i);
	if i<4
		s2=s2+i*b(i)*sqrt((M1+i*M2)/(i*(M1+M2))) ;
	end
end
y=s2/s1;
	
