# first passage set up to compute the firing times
# -1 = f(V) T'(V) + sigma^2/2 T''(V)   -oo < V < Vspike
# we desire T(Vreset) but we cannot get this directly
# from AUTO, so we will instead define two variables u,w 
# defined on [0,1] and scaled so that u(1)=w(0) = T(Vreset)
# 
# 
# this is a bit weird, but here goes:
# u is defined on -a < V < Vreset
# a is sort of like infinity. b=Vreset+a is the length of this interval
# w is defined on  Vreset < V < VSpike  c = (Vspike-Vreset) is length of this
# u(s=1)=w(s=0)
# note d/dV = d/d(bs) for b=(Vreset+a) 
#      d/dV = d/d(cs) for c=(Vspike-Vreset)
# c u'(s=1)=b w'(s=0)
# now what about the BC at V=-a.  
# easiest is to make T'(-a)=0
# or   u'(s=0)=0
#
# for better accuracy, T'(V) f(V) ~ -1 for V large
# so for QIF T'(V)= -1/a^2
# and thus   u'(s=0) = -b/a^2
#
# clearly
# w(s=1)=0 since T(Vspike)=0
# 
# this is set up for sigma=1, I=50 easy to solve! 
# for I large and negative or sigma small, XPP cannot solve it
# using shooting, but AUTO can. Below it is set up for AUTO
# with struggle AUTO will take it down to almost  -2
#
# in XPP:
# Click on Bndry Val Show
# File AUTO
# in AUTO 
# Run BoundryValue
# youll get the FI curve
# QIF:
f(x)=x^2+I
init u=1,w=1,wp=-2,freq=1 
par I=50,sig=1,vreset=-10,vspike=10,a=20
!b=(vreset+a)
!c=(vspike-vreset)

du/dt=up
dup/dt=-b*b/sig-f(-a+b*s)*up*b/sig
dw/dt=wp
dwp/dt=-c*c/sig-f(vreset+c*s)*wp*c/sig
ds/dt=1
dfreq/dt=0
b up
b w'
b c*up'-b*wp
b u'-w
b s
b freq-1/w
# I add the actual values of V here to plot the whole thing
# on the real V scale
aux vleft=-a+b*s
aux vright=vreset+c*s
# numerics stuff
@ total=1,dt=.005,meth=cvode,tol=1e-8,atol=1e-8
# plot both segments
@ nplot=2
@ yp=u,yp2=w,xp=vleft,xp2=vright,xlo=-10,xhi=10,ylo=0,yhi=10
# AUTO parameters
@ NTST=60,DSMIN=1e-4,NPR=200,NMAX=1000,DS=-.01
@ PARMIN=-2,PARMAX=60
# Assumed here that I will be the parameter to get "f-I" curve
# freq=1/w(0) will be plotted - this is 1/T(Vreset)
@ AUTOXMIN=-2,AUTOXMAX=50,AUTOYMAX=10,AUTOVAR=freq,AUTOYMIN=0
done