wf = function(npop=5,ngen=500,nall=100,init=1,mu=0)
# this is a comment. Comments inform users; they are ignored by R
# this function simulates the frequency of one allele with initial frequency init in each of npop haploid populations each of size nall, which then evolves independently in each population over ngen generations according to a haploid Wright-Fisher model with mutation rate mu per generation
{
	freq = matrix(init,npop,ngen) # create a matrix with npop rows and ngen columns, initialise every entry to init
	for(i in 1:(ngen-1)) freq[,i+1] = rbinom(npop,nall,mu+(1-2*mu)*freq[,i]/nall) # create the 2,3,4, ... ,ngen generations, in each of the npop populations, by random sampling from the binomial distribution.  The expected allele fraction in the new generation is the allele fraction in the current generation (freq[,i]/nall) modified by the effect of mutation

# by commenting out one of the following two lines, you can choose only to plot the frequency in each generation (1 plot with a different line for each population), or only to output the final-generation frequencies in a vector
 	matplot(t(freq),type="l",lty=1,lwd=2,xlab="generations",ylab="allele count",main=paste("WF model in", npop," populations, mu = ",mu))
# 	return(freq[,ngen])
}