function p = OccupationProbability(epsilon,K)

global pE pW pN pS;
pE = 0.25+epsilon; pN = 0.25; 
pW = 0.25-epsilon; pS = 0.25;

p = zeros(K,1);
p(1) = 1; Pi = 1; 
for k=1:K-1
    Pi = step(step(Pi));
    p(k+1) = Pi(k+1,k+1);
end
    
return;

function PiNew = step(Pi)

global pE pW pN pS;

[k,k] = size(Pi); PiNew = zeros(k+1,k+1); i=1:k;

PiNew(i+1,i+1) = pE*Pi;
PiNew(i  ,i  ) = PiNew(i  ,i  ) + pW*Pi;
PiNew(i  ,i+1) = PiNew(i  ,i+1) + pN*Pi;
PiNew(i+1,i  ) = PiNew(i+1,i  ) + pS*Pi;

return;