P=[]
G=10;
N=80;
%N=1;
Q=[];
for j=1:N
a=rand(1,3);
%a=[0 1 0] % for Mendel's 3:1 ratio rule. 
a(1)=(j-1)/(N);
a(2:3)=a(2:3)*(1-a(1))/(a(2)+a(3));
%p=a./sum(a);
p=a;
P=[p];
for i=1:G
    q=peaplants(p);
    q=q./sum(q);
    Q=[Q;q];
    P=[P;q];
    p=q;
end
[P(1,:); P(end,:)]
end

figure(1)
hold off

plot(Q(:,1),Q(:,2),'ko',1-Q(:,1)-Q(:,2),Q(:,2),'ro')
xlabel('P_{YY}')
ylabel('P_{Yg}')
title('Steady State Distribution')

figure(2)
hold off
plot3(Q(:,1),Q(:,2),1-Q(:,1)-Q(:,2),'ko',1-Q(:,1)-Q(:,2),Q(:,2),Q(:,1),'ro')
xlabel('P_{YY}')
ylabel('P_{Yg}')
zlabel('P_{gg}')
title('Steady State Distribution')
grid on

return
%%
a=rand(1,3);
p=a./sum(a);
q=peaplants(p);
q
sum(q)

%%