function output = EhrenfestPower(beads , steps )
%EhrenfestPower exponentiates the basic Ehrenfest model matrix

EhrenfestMatrix=zeros(beads+1);

for i=1:beads+1
    if i==1
        EhrenfestMatrix(i,i+1)=1;
    elseif i==beads+1
        EhrenfestMatrix(i,i-1)=1;
    else
        EhrenfestMatrix(i,i-1)=(i-1)/beads;
        EhrenfestMatrix(i,i+1)=(beads-i+1)/beads;
    end

output=EhrenfestMatrix^steps;

end

