# lattice prisoners dilemma
@ meth=discrete,total=1001
par dt=.2
par rr=3,ss=0,tt=5,pp=1
par w=.6
par alpha=.2
par z=2
!vtt=rr/(1-w)
!vtd=ss+w*pp/(1-w)
!vdt=tt+w*pp/(1-w)
!vdd=pp/(1-w)
mt(n)=exp(-alpha*(1-w)*(n*vtt+(z-n)*vtd))
md(n)=exp(-alpha*(1-w)*(n*vdt+(z-n)*vdd))
x0=x500
x501=x1
# n=0 or 1 - left or right
n[1..500]=ran(1)<.5
# dth=1 if you die and 0 if you live
dth[1..500]=heav((x[j]*mt(x[j+1]+x[j-1])+(1-x[j])*md(x[j+1]+x[j-1]))*dt-ran(1))
x[1..500]'=x[j]*(1-dth[j])+dth[j]*(n[j]*x[j-1]+(1-n[j])*x[j+1])
x[1..500](0)=ran(1)<.5
aux rhot=sum(0,499)of(shift(x1,i'))/500
done