Imports System.Math
Imports System.IO
Module Module1
 Sub Main()
  Dim auc As Double = 0
  Dim ba As Double = 1
  Dim cl As Double = 0
  Dim cs1, cs2, cs3, cs4, cs5, cs6, cs7, cs8 As Double
  cs1 = 0.5 : cs2 = 0.5 : cs3 = 0.5 : cs4 = 0.5 : cs5 = 0.5 : cs6 = 0.5 : cs7 = 0.5 : cs8 = 0.5
  Dim d As Double = 0.0003
  Dim dlt As Double = 0.003
  Dim dose As Double = 0.680678408
  Dim f(4, 8) As Double
  Dim g As Double = 0
  Dim h As Double = 0.0001
  Dim i, j As Integer
  Dim ka1, ka2, ka3, ka4, ka5, ka6, ka7 As Double
  ka1 = 0.03 : ka2 = 0.03 : ka3 = 0.03 : ka4 = 0.03 : ka5 = 0.03 : ka6 = 0.03 : ka7 = 0.03
  Dim kn1(4, 8), kn2(4, 8), kn3(4, 8), kn4(4, 8) As Double
  Dim k12 As Double = 0
  Dim k21 As Double = 0
  Dim kt1, kt2, kt3, kt4, kt5, kt6, kt7 As Double
  kt1 = 0.0352 : kt2 = 0.0352 : kt3 = 0.0352 : kt4 = 0.0352 : kt5 = 0.0352 : kt6 = 0.0352 : kt7 = 0.0352
  Dim p As Double = 1300
  Dim pc As Double
  Dim q As Double = 5
  Dim r As Double = 0.05
  Dim s(4) As Double
  Dim t As Double = 0.05
  Dim V1, V2, V3, V4, V5, V6, V7, V8 As Double
  V1 = 35.7 : V2 = 35.7 : V3 = 35.7 : V4 = 35.7 : V5 = 35.7 : V6 = 35.7 : V7 = 35.7 : V8 = 35.7
  Dim vd As Double = 1000
  Dim w As Double = 70
  Dim X As Double = 0
  Dim xmax As Double = 720
  Dim Xso As Double = 10
  Dim Y(4, 8) As Double
  Dim yn(4, 8) As Double
  Dim file As StreamWriter
  file = My.Computer.FileSystem.OpenTextFileWriter("C:\Users\Kevin C. Johnson\cat.txt", False)
  For i = 1 To 4
   For j = 1 To 8
    Y(i, j) = 0
   Next
  Next
  Y(1, 1) = dose
  Dim col As Double
Start:
  col = Y(1, 8) / dose * 100
  If X > 0.1 Then h = 0.001
  If X > 1 Then h = 0.01
  For i = 1 To 4
   s(i) = 0
  Next
  For i = 1 To 4
   For j = 1 To 8
    s(i) = s(i) + Y(i, j)
   Next
  Next
  t = r
  If t > dlt Then t = dlt
  pc = s(2) / dose * 100
  auc = auc + pc * h
  If g < h Then file.Write(X.ToString("0.00"))
  If g < h Then file.Write(" ")
  If g < h Then file.WriteLine(s(3) / dose * 100.ToString("0.00E+0"))
  g = g + h
  If g > q - h / 2 Then g = 0
  For i = 1 To 4
   For j = 1 To 8
    eq(f, d, Y, p, t, r, cs1, cs2, cs3, cs4, cs5, cs6, cs7, cs8,
       V1, V2, V3, V4, V5, V6, V7, V8, kt1, kt2, kt3, kt4, kt5, kt6, kt7,
       ka1, ka2, ka3, ka4, ka5, ka6, ka7, cl, vd, k12, k21, ba)
    yn(i, j) = Y(i, j)
    kn1(i, j) = f(i, j)
    X = X + h / 2
    Y(i, j) = yn(i, j) + kn1(i, j) * h / 2
    If Y(1, j) < 1.0E-17 Then Y(1, j) = 0
    eq(f, d, Y, p, t, r, cs1, cs2, cs3, cs4, cs5, cs6, cs7, cs8,
       V1, V2, V3, V4, V5, V6, V7, V8, kt1, kt2, kt3, kt4, kt5, kt6, kt7,
       ka1, ka2, ka3, ka4, ka5, ka6, ka7, cl, vd, k12, k21, ba)
    kn2(i, j) = f(i, j)
    Y(i, j) = yn(i, j) + kn2(i, j) * h / 2
    If Y(1, j) < 1.0E-17 Then Y(1, j) = 0
    eq(f, d, Y, p, t, r, cs1, cs2, cs3, cs4, cs5, cs6, cs7, cs8,
       V1, V2, V3, V4, V5, V6, V7, V8, kt1, kt2, kt3, kt4, kt5, kt6, kt7,
       ka1, ka2, ka3, ka4, ka5, ka6, ka7, cl, vd, k12, k21, ba)
    kn3(i, j) = f(i, j)
    X = X + h / 2
    Y(i, j) = yn(i, j) + kn3(i, j) * h
    If Y(1, j) < 1.0E-17 Then Y(1, j) = 0
    eq(f, d, Y, p, t, r, cs1, cs2, cs3, cs4, cs5, cs6, cs7, cs8,
       V1, V2, V3, V4, V5, V6, V7, V8, kt1, kt2, kt3, kt4, kt5, kt6, kt7,
       ka1, ka2, ka3, ka4, ka5, ka6, ka7, cl, vd, k12, k21, ba)
    kn4(i, j) = f(i, j)
    Y(i, j) = yn(i, j) + h * (kn1(i, j) + 2 * kn2(i, j) + 2 * kn3(i, j) + kn4(i, j)) / 6
    If Y(1, j) < 1.0E-17 Then Y(1, j) = 0
    X = X - h
   Next
  Next
  X = X + h
  If X < xmax + h Then GoTo Start
  file.Close()
 End Sub
 Sub eq(ByRef f(,) As Double, ByRef d As Double, ByRef Y(,) As Double, ByRef p As Double, ByRef t As Double,
        ByRef r As Double, ByRef cs1 As Double, ByRef cs2 As Double, ByRef cs3 As Double, ByRef cs4 As Double,
        ByRef cs5 As Double, ByRef cs6 As Double, ByRef cs7 As Double, ByRef cs8 As Double, ByRef V1 As Double,
        ByRef V2 As Double, ByRef V3 As Double, ByRef V4 As Double, ByRef V5 As Double, ByRef V6 As Double,
        ByRef V7 As Double, ByRef V8 As Double, ByRef kt1 As Double, ByRef kt2 As Double, ByRef kt3 As Double,
        ByRef kt4 As Double, ByRef kt5 As Double, ByRef kt6 As Double, ByRef kt7 As Double, ByRef ka1 As Double,
        ByRef ka2 As Double, ByRef ka3 As Double, ByRef ka4 As Double, ByRef ka5 As Double, ByRef ka6 As Double,
        ByRef ka7 As Double, ByRef cl As Double, ByRef vd As Double, ByRef k12 As Double, ByRef k21 As Double,
        ByRef ba As Double)
  f(1, 1) = -3 * d * Y(1, 1) / p / t / r * (cs1 - Y(2, 1) / V1) - kt1 * Y(1, 1)
  f(2, 1) = 3 * d * Y(1, 1) / p / t / r * (cs1 - Y(2, 1) / V1) - (kt1 + ka1) * Y(2, 1)
  f(3, 1) = ka1 * Y(2, 1) * ba - (cl / vd + k12) * Y(3, 1) + k21 * Y(4, 1)
  f(4, 1) = k12 * Y(3, 1) - k21 * Y(4, 1)
  f(1, 2) = kt1 * Y(1, 1) - 3 * d * Y(1, 2) / p / t / r * (cs2 - Y(2, 2) / V2) - kt2 * Y(1, 2)
  f(2, 2) = kt1 * Y(2, 1) + 3 * d * Y(1, 2) / p / t / r * (cs2 - Y(2, 2) / V2) - (kt2 + ka2) * Y(2, 2)
  f(3, 2) = ka2 * Y(2, 2) * ba - (cl / vd + k12) * Y(3, 2) + k21 * Y(4, 2)
  f(4, 2) = k12 * Y(3, 2) - k21 * Y(4, 2)
  f(1, 3) = kt2 * Y(1, 2) - 3 * d * Y(1, 3) / p / t / r * (cs3 - Y(2, 3) / V3) - kt3 * Y(1, 3)
  f(2, 3) = kt2 * Y(2, 2) + 3 * d * Y(1, 3) / p / t / r * (cs3 - Y(2, 3) / V3) - (kt3 + ka3) * Y(2, 3)
  f(3, 3) = ka3 * Y(2, 3) * ba - (cl / vd + k12) * Y(3, 3) + k21 * Y(4, 3)
  f(4, 3) = k12 * Y(3, 3) - k21 * Y(4, 3)
  f(1, 4) = kt3 * Y(1, 3) - 3 * d * Y(1, 4) / p / t / r * (cs4 - Y(2, 4) / V4) - kt4 * Y(1, 4)
  f(2, 4) = kt3 * Y(2, 3) + 3 * d * Y(1, 4) / p / t / r * (cs4 - Y(2, 4) / V4) - (kt4 + ka4) * Y(2, 4)
  f(3, 4) = ka4 * Y(2, 4) * ba - (cl / vd + k12) * Y(3, 4) + k21 * Y(4, 4)
  f(4, 4) = k12 * Y(3, 4) - k21 * Y(4, 4)
  f(1, 5) = kt4 * Y(1, 4) - 3 * d * Y(1, 5) / p / t / r * (cs5 - Y(2, 5) / V5) - kt5 * Y(1, 5)
  f(2, 5) = kt4 * Y(2, 4) + 3 * d * Y(1, 5) / p / t / r * (cs5 - Y(2, 5) / V5) - (kt5 + ka5) * Y(2, 5)
  f(3, 5) = ka5 * Y(2, 5) * ba - (cl / vd + k12) * Y(3, 5) + k21 * Y(4, 5)
  f(4, 5) = k12 * Y(3, 5) - k21 * Y(4, 5)
  f(1, 6) = kt5 * Y(1, 5) - 3 * d * Y(1, 6) / p / t / r * (cs6 - Y(2, 6) / V6) - kt6 * Y(1, 6)
  f(2, 6) = kt5 * Y(2, 5) + 3 * d * Y(1, 6) / p / t / r * (cs6 - Y(2, 6) / V6) - (kt6 + ka6) * Y(2, 6)
  f(3, 6) = ka6 * Y(2, 6) * ba - (cl / vd + k12) * Y(3, 6) + k21 * Y(4, 6)
  f(4, 6) = k12 * Y(3, 6) - k21 * Y(4, 6)
  f(1, 7) = kt6 * Y(1, 6) - 3 * d * Y(1, 7) / p / t / r * (cs7 - Y(2, 7) / V7) - kt7 * Y(1, 7)
  f(2, 7) = kt6 * Y(2, 6) + 3 * d * Y(1, 7) / p / t / r * (cs7 - Y(2, 7) / V7) - (kt7 + ka7) * Y(2, 7)
  f(3, 7) = ka7 * Y(2, 7) * ba - (cl / vd + k12) * Y(3, 7) + k21 * Y(4, 7)
  f(4, 7) = k12 * Y(3, 7) - k21 * Y(4, 7)
  f(1, 8) = kt7 * Y(1, 7) - 3 * d * Y(1, 8) / p / t / r * (cs8 - Y(2, 8) / V8)
  f(2, 8) = kt7 * Y(2, 7) + 3 * d * Y(1, 8) / p / t / r * (cs8 - Y(2, 8) / V8)
 End Sub
End Module