N <- 2000;
num_trials <- 1000;
net_profits <- list();
for(tn in 1:num_trials) {
    coinTossResults   <- cumsum(sample(c(-1,1), N, replace = TRUE))
    net_profits[[tn]] <- coinTossResults[length(coinTossResults)]
}

x <- unlist(net_profits)
h<-hist(x, breaks=40, col="red", xlab="Profit",
   main="Histogram of profit/loss distribution")
xfit<-seq(min(x),max(x),length=length(x))
yfit<-dnorm(xfit,mean=mean(x),sd=sd(x))
yfit <- yfit*diff(h$mids[1:2])*length(x)
lines(xfit, yfit, col="blue", lwd=2)