% simulate two tank system

function [yt,theta]=system(y,u,t,my)

m=2;
% transition probability matrix-theta

for j=1:my^2
   theta(1,j)=(1-j/(my^m+1))^(my-1);
   for i=2:my
       i1=i-1;
       deno=1;
       for k=my-i1:my-1
         deno=deno*k;
       end
       numo=1;
       for k=1:i1
         numo=numo*k;
       end
       comb=deno/numo;
       theta(i,j)=comb*(j/(my^m+1))^i1*(1-j/(my^m+1))^(my-1-i1);
   end
end   

% generating data y(t)
ty=0;                    

  ty=ty+(y(t-1)-1)+my*(u(t)-1);

ty=ty+1;


  rn=rand;
  pom=0;
  for f=1:my
     if (pom<rn)&(pom+theta(f,ty)>rn)
       yt=f;               
     end
     pom=pom+theta(f,ty);
  end