! Ten, version 2.
! John Wise, initiated 2001/01/12, last updated 2001/01/18.
! Based on "Decay" in "Computational Physics" by N. Giordano.
program Ten
   library "IrGraphics*"
   option nolet
   dim t(0), na(0), nb(0), nc(0), nd(0)
   call Init(nt, dt, t0, na0, nb0, nc0, nd0, ra, rb)
   mat redim t(0:nt), na(0:nt), nb(0:nt), nc(0:nt), nd(0:nt)
   t(0) = t0
   na(0) = na0
   nb(0) = nb0
   nc(0) = nc0
   nd(0) = nd0
   call Calc(nt, dt, t(), na(), nb(), nc(), nd(), ra, rb)
   id$ = "JHW: Problem 10. X-c : Dashed Line. X-d : Solid Line"
   call Graph2A(id$, "Time", "Number", t(), nc(), 2, t(), nd(), 1)
   print " t = "; t(nt)
   print " nc = "; nc(nt)
   print " nd = "; nd(nt)
   print "Done"
end
sub Init(nt, dt, t0, na0, nb0, nc0, nd0, ra, rb)
   nt = 100           ! number of time steps
   dt = 0.01          ! time step
   t0 = 0             ! initial time
   na0 = 100           ! initial value
   nb0 = 10             ! initial value
   nc0 = na0*nb0       ! initial value
   ra = 5               ! decay rate (a)
   rb = 2               ! decay rate (b)
end sub
sub Calc(nt, dt, t(), na(), nb(), nc(), nd(), ra, rb)
   for i = 1 to nt
      nd1 = -ra*na(i-1)
      nd2 = -rb/nb(i-1)
      nd3 = -(ra+rb)*nc(i-1)
      t(i) = t(i-1) + dt
      na(i) = na(i-1) + nd1*dt
      nb(i) = nb(i-1) + nd2*dt
      nc(i) = nc(i-1) + nd3*dt
      nd(i) = na(i)*nb(i)
   next i
end sub