% urn model
iter=500;
real=1;
store=zeros(real,iter+1);

for k=1:real

    
% condizione iniziale.

A=[2,0];
B=[0,3];
store(k,1)=1;


for i=1:iter
   
    p1=rand(1);
    p2=rand(1);
    pA=A(1)/2;
    pB=B(1)/3;
    
    if p1<pA
        A(1)=A(1)-1;
        B(1)=B(1)+1;
    else
        A(2)=A(2)-1;
        B(2)=B(2)+1;
    end
    
    if p2<pB
        A(1)=A(1)+1;
        B(1)=B(1)-1;
    else
        A(2)=A(2)+1;
        B(2)=B(2)-1;
    end
    
    if A(1)==2
    store(k,i+1)= 1;
    elseif A(1)==1
    store(k,i+1)=2;
    elseif A(1)==0
     store(k,i+1)=3;
    end 
figure(1)
subplot(2,1,1)
plot(1,A(1),'ko','MarkerSize',20)
hold on
plot(2,A(2),'ro','MarkerSize',20)
plot([0.5 2.5],[1 1],'k--')
plot([0.5 2.5],[2 2],'k--')
hold off
axis([0.5 2.5 0 2.5])
subplot(2,1,2)
plot(1,B(1),'ko','MarkerSize',20)
hold on
plot(2,B(2),'ro','MarkerSize',20)
plot([0.5 2.5],[1 1],'k--')
plot([0.5 2.5],[2 2],'k--')
plot([0.5 2.5],[3 3],'k--')
hold off
axis([0.5 2.5 0 2.5])
pause(0.1)


end

figure(2)
[Y X]=hist(store(k,:),[1 2 3]);
bar(X,Y/length(store(k,:)))
hold on 
plot([0.5 3.5],[1/10 1/10],'k--')
plot([0.5 3.5],[3/5 3/5],'k--')
plot([0.5 3.5],[3/10 3/10],'k--')
hold on
end



figure(3)
plot(mean(store),'-o')
%plot(mean(store),'-o') %if real=1
hold on 
for i=0:iter
plot(i+1,11/5-6/5*(1/6)^i,'r*')
end