      PROGRAM plot_function                 ! plot analytical solution
      IMPLICIT NONE
      REAL*8 T0, Ts, r, t, tmax
      INTEGER IR_
      CALL GWopen(IR_, 0)
      CALL initial(t,tmax,r,Ts,T0)
      CALL set_up_window(REAL(t),REAL(tmax),REAL(Ts),REAL(T0))
      CALL show_plot(t,tmax,r,Ts,T0)
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(t,tmax,r,Ts,T0)
      IMPLICIT NONE
      REAL*8 t, tmax, r, Ts, T0
      t = 0
      T0 = 83                               ! initial coffee temperature (C)
      Ts = 22                               ! room temperature (C)
      WRITE(*,'(A,$)') 'cooling constant r = '
      READ(*,*) r                           ! cooling constant
      WRITE(*,'(A,$)') 'duration = '
      READ(*,*) tmax                        ! time (minutes)
      END
C
      SUBROUTINE set_up_window(xmin,xmax,ymin,ymax)
      IMPLICIT NONE
      REAL mx, my, xmin, xmax, ymin, ymax
      INTEGER IR_
      mx = 0.01*(xmax - xmin)               ! margin
      my = 0.01*(ymax - ymin)
      CALL GWindow(IR_, xmin - mx, ymin - my, xmax + mx, ymax + my)
*     default background color black on all computers except Macintosh
      CALL GWmove2(IR_, xmin, ymin)         ! abbreviation for PLOT LINES:
      CALL GWline2(IR_, xmax, ymin)
      CALL GWline2(IR_, xmin, ymin)
      CALL GWline2(IR_, xmin, ymax)
      END
C
      SUBROUTINE show_plot(t,tmax,r,Ts,T0)
      IMPLICIT NONE
      REAL*8 t, tmax, r, Ts, T0, f
      INTEGER IR_
      CALL GWsetpen(IR_, -1, -1, -1, -1)
      CALL GWsetpen(IR_, 13, -1, -1, -1)
      DO WHILE(t .LE. tmax)
         CALL GWline2(IR_, REAL(t), REAL(f(t,r,Ts,T0))) ! abbreviation for PLOT LINES
         t = t + 0.1D0
      END DO
      END
C
      FUNCTION f(t,r,Ts,T0)
      IMPLICIT NONE
      REAL*8 f, t, r, Ts, T0
      f = Ts - (Ts - T0)*exp(-r*t)
      END
