# d'Onofrio 2005

# parametres :
p alpha=1

# Kuznetsov et al.
f(x) = 1.636*(1 - 0.002*x)/alpha
phi(x) = 1
beta(x) = (1.131*x)/(20.19 + x)

# sigma*q(x)=0.1181 -> q(x)=0.1181/sigma ?
# sigma n'apparait pas explicitement, on l'absorbe dans s_q
s_q(x) = 0.1181

# tumeur non-aggressive :
#mu(x) = 0.00311*x + 0.3743

#tumeur agressive :
mu(x) =10*(0.00311*x) + 0.3743

# equation (8)
dx/dt = x*(alpha*f(x) - phi(x)*y)
# equation (9)
dy/dt = -(mu(x) -beta(x))*y + s_q(x) + theta(t)

# traitement

p A=0.50
p b=0,omega=0
theta(t)=A*(1+b*cos(omega*t))

init x=0.1,y=0.6

# condition de croissance
aux xp=if(x>100)then(omega)else(-100)
aux yp=if(x>100)then(b)else(-100)



@ meth=qualrk,total=1000 #total = temps
@ maxstor=1000000
@ bound=1000

# axes figure: 
@ xp=xp, yp=yp, xlo=0,ylo=0, xhi=0.5, yhi=1,lt=-1


# AUTO bifurcation settings
# Parametre de bifurcation : theta_m 
# Diagramme : X vs theta_m
@autoxmin=0,autoxmax=2,autoymin=0, autoymax=600

@dsmax=0.7,dsmin=.001, parmin=-2,parmax=2.0, ds=0.02