# Bistable wave simulation

# time in s
# space in um
# d in um^2/s
# parameters
param d=30,dx=0.2,sigma=5,taur=1,taud=10000000,cth=0.1 

#the initial condition, equation and a function
c[0..100](0)=0

s[0..20](0)=1.0
s[21..100](0)=0

#global 1 c0-cth {s0=1}
global 1 c10-cth {s10=1}
global 1 c20-cth {s20=1}
global 1 c30-cth {s30=1}
global 1 c40-cth {s40=1}
global 1 c50-cth {s50=1}
global 1 c60-cth {s60=1}
global 1 c70-cth {s70=1}
global 1 c80-cth {s80=1}
global 1 c90-cth {s90=1}
global 1 c100-cth {s100=1}

dc0/dt = sigma*heav(s0)/2/taur-c0/taud+d*(c1-c0)/dx^2
dc[1..99]/dt = sigma*heav(s[j])/taur-c[j]/taud+d*(c[j-1]-2*c[j]+c[j+1])/dx^2
dc100/dt = sigma*heav(s100)/2/taur-c100/taud+d*(c99-c100)/dx^2

ds[0..100]/dt = -heav(s[j])/taur

aux logc[0..100] = c[j] 

#program end
@ total=2,trans=0,dt=0.0001,xlo=0,xhi=2000,ylo=0,yhi=1
@ maxstore=1000000,bounds=10000
@ xplot=x,yplot=Ca10
done