@ autoeval=0
# generate Lorenzian of frequencies
table w % 400 0 399  tan(pi*(ran(1)-.5))
par tau=.5,vspike=100,mu=5,del=.25,g=0
par dt=.005
# synapse
#  fr is fraction that jumped past threshold
s'=s -s*dt/tau+fr/tau
# real time
tt'=tt+dt
# update vnew
vnew[0..399]=v[j]+dt*(v[j]^2+mu+del*w([j])+g*s)
# is vnew > vspike?
flag[0..399]=vnew[j]>vspike
# count fraction of guys that have jumped
fr=sum(0,399)of(shift(flag0,i'))/400
# update v with appropriate reset
v[0..399]'=if(flag[j])then(-vspike)else(vnew[j])
# initial conditions
v[0..399](0)=-50+ran(1)*100
# average potential
aux vbar=sum(0,399)of(shift(v0,i'))/400
# keep track of average firing rate
aux fr=fr
@ bound=1000
# nice plots
@ xp=tt,yp=s,xlo=0,xhi=200,ylo=0,yhi=1.25
@ total=40000,meth=discrete,nout=2
d