! Population, version 1.
! John Wise, initiated 2001/01/12, last updated 2001/01/18.
! Based on "Decay" in "Computational Physics" by N. Giordano.
program Population
   library "IrGraphics*"
   option nolet
   id$ = "JHW: Population, n0:100, Birth:1/20, Death:3/10000"
   dim t(0), n(0)
   call Init(nt, dt, t0, n0, a, b)
   mat redim t(0:nt), n(0:nt)
   t(0) = t0
   n(0) = n0
   call Calc(nt, dt, t(), n(), a, b)
   call Graph1A(id$, "Time", "Population", t(), n(), 1)
   print " t = "; t(nt)
   print " n = "; n(nt)
   print "Done"
end
sub Init(nt, dt, t0, n0, a, b)
   nt = 200          ! number of time steps
   dt = 0.5          ! time step
   t0 = 0            ! initial time
   n0 = 100            ! initial value
   a = 5/n0            ! birth rate
   b = 3/(n0*n0)       ! death rate
end sub
sub Calc(nt, dt, t(), n(), a, b)
   for i = 1 to nt
      nd = a*n(i-1) - b*n(i-1)*n(i-1)
      t(i) = t(i-1) + dt
      n(i) = n(i-1) + nd*dt
   next i
end sub