share/contrib/Eulix/Eulix_T5N.mac

   1 load("Eulix.mac")$
   2
   3 Start:elapsed_real_time()$
   4
   5 n:1000$ p:10$
   6 rx:2.4E-5$ rv:7.4E-4$
   7 sx:0.02$ sy:0.265$
   8 ginf:0.01$
   9 vg:100$ sv:0$
  10 g(v):= ginf*v/(v+100)$
  11 f(v):= 20.0*v/(v+1.0)$
  12
  13 debugmode(true)$
  14 /* :break Eulix$ */
  15 t_final: 12*365+4$
  16
  17 Logging:false$
  18
  19 [Root_at,y_root]: Eulix([p - sx*x -rx*x*v +g(v)*x,
  20                           rx*x*v -sy*y -g(v)*y,
  21                           n*g(v)*y +f(v) -sv*v -(rx+rv)*x*v],[x,y,v],
  22                           [1000.0,0,1.0],[t,0.0,t_final,1,v-t-245.0],
  23                           initial_step_size=1,absolute_tolerance=1e-8,logging=Logging,
  24                           'tabular='none)$
  25 if first(Root_at) > 0 then (
  26   print("I've used ",elapsed_real_time()-Start,"seconds for computing the AIDS model"),
  27   print("Solution triggered by condition",first(Root_at)," at time ",second(Root_at),
  28         "  and state:",third(Root_at)),
  29   print("y_root=",y_root),
  30   t1:second(Root_at),
  31   y_final: Eulix([p - sx*x -rx*x*v +g(v)*x,
  32                   rx*x*v -sy*y -g(v)*y,
  33                   n*g(v)*y +f(v) -sv*v -(rx+rv)*x*v],
  34                   [x,y,v],y_root,[t,t1,t_final,1],
  35                   absolute_tolerance=1e-8,logging=true,
  36                   'tabular='none),
  37
  38   print("solution at final time",t_final," is ",y_final)
  39 ) else print("no root condition has been triggered")$
  40