# Bistable wave simulation

#parameters
param ip3=.7, a2=.2, Caer=1
param d1=0.1, d2=2, d3=0.2, d5=0.2
param v1=20, v2=0.004, v3=1.2, k3=0.15
param deff=16,dx=10
# dx in units of um
# deff in units of um^2/sec
#the initial condition, equation and a function
Ca[0..4](0)=.2
Ca[5..100](0)=.01

dCa0/dt = f(Ca0)+deff*(Ca1-Ca0)/dx^2
dCa[1..99]/dt = f(Ca[j])+deff*(Ca[j-1]-2*Ca[j]+Ca[j+1])/dx^2
dCa100/dt = f(Ca100)+deff*(Ca99-Ca100)/dx^2

f(Ca) = (v2+v1*(winf(Ca)*Ca/(Ca+d5))^3)*(Caer-Ca)-v3*Ca^2/(k3^2+Ca^2)
winf(Ca) = (d2*ip3/(ip3+d3))/(Ca+(d2*(ip3+d1)/(ip3+d3)))
 
#program end
@ total=200,trans=0,DT=.1,xlo=0,xhi=200,ylo=0,yhi=1
@ maxstore=1000000,bounds=10000
@ xplot=x,yplot=Ca0
done