function x = runepidemic(p,q,N,T)
    x=zeros(T,1); %Number of individuals infected initially.
    x(1)=1;

    %For all time steps
    for t=1:T
        %Find the number recovering per time step using binomial distribution
        recover=binornd(x(t),p);
        %Find the number infecting others using binomial distribution
        contacted=binornd(x(t),q);
        
        %Now find number of contacted indivduals who are not infected already.
        %Here we have to be careful to take in to account of the fact that
        %if the same individual is contacted twice, they are only infected
        %once. We first make a random list of all the individuals contacted
        contacts=ceil(rand(contacted,1)*N);
        %Now find the unique set of contacted individuals
        contacts_unique=unique(contacts);
        %Then finally consider only those contacts which are not infected 
        %already 
        contacted_s=sum(contacts_unique>x(t));
        %Update infected population
        x(t+1)=x(t)-recover+contacted_s;
    end