function [] = PlotSimulate(A, varargin)

% The next line sets up Matlab's random number generator.
rand('twister', sum(100*clock)); 

n = size(A, 1);
Trials = 1000; % Trials is the number of random users.
Steps = 50; % Steps is the number of clicks our random users will make.

Count = zeros(n, Steps+1); % Count will be used to store the number of users 
                          % on each page at each time.
Deg = sum(A); % Deg stores the number of links from each page.

B=A; % We make a copy of B to preserve A, in case we wish to use it 
     % for some other purpose.
     
     
% The next lines attempt to turn each column of the adjacency matrix into a
% probability vector, by dividing by the number of links leaving the 
% corresponding vertex. If a page does not link to any other page this would 
% lead to dividing by zero, so whenever this happens we keep the column of
% zeros.

for i=1:n
    if (Deg(i) > 0)
        B(:, i)=B(:, i)./Deg(i);
    else 
        B(:, i)=0;
    end
end

% The next lines perform the simulation.

for i=1:Trials
    Position = randi(n); % We start each surfer at a randomly chosen vertex.
                         % We use Position to record a surfer's location.
    Count(Position, 1) = Count(Position, 1)+1; 
    for j=2:Steps+1
        % If there is a link leaving a verex then we update Position to
        % reflect a random choice of link to follow. Otherwise we do not
        % change Position and the surfer is stuck on the page where they are.
        if Deg(Position) > 0 
            Position = ProbGen(B(:, Position));
        end
    Count(Position, j) = Count(Position, j)+1;    
        % Each time we update the Position of a surfer we add one to the
        % Count of users in that Position at that time.
    end
end
Probabilities = Count ./ Trials
        % To turn Count into a set of Probability distributions, we divide
        % every element in the Count matrix by Trials.

% The next lines plot how the probabilities of being in certain positions
% specified by the surfer evolve over time. 
        
nArgs = length(varargin);
hold on

axis([1 Steps+1 0 1])

for i=1:nArgs
    plot (Probabilities(varargin{i}, :))
end
hold off
end





