share/contrib/Eulix/Eulix_TS.mac

   1 load("Eulix.mac")$
   2 load("rkf45.mac")$
   3
   4 /* This is a 1-d example
   5    Shampine's ball of flame
   6    y'= y^2-y^3   y(0)= delta   [0,1/delta]
   7 */
   8
   9 debugmode:true$
  10 delta:0.00001$
  11 timer(rkf45,Eulix)$
  12
  13 Logging:false$
  14
  15 RR:rkf45(y^2-y^3,y,delta,[t,0,2/delta],
  16          report=true,max_iterations=40000,absolute_tolerance=5E-8)$
  17
  18 [tlist,ylist]: Eulix(y^2-y^3,y,delta,[t,0,2/delta],
  19                      logging=Logging,absolute_tolerance=5E-8,combined_t_y_list=false)$
  20 timer_info();  /* rkf45: 37.376 sec   Eulix: 2.192 sec */
  21
  22 plot2d([discrete,tlist,map(lambda([z],z[1]),ylist)],
  23        [legend,"y"],[title,"Shampine's ball of flame"])$
  24 read("enter ; or $ to quit");
  25
  26