VERSION 5.00
Begin VB.Form Form2 
   Caption         =   "Form2"
   ClientHeight    =   8595
   ClientLeft      =   60
   ClientTop       =   120
   ClientWidth     =   11880
   LinkTopic       =   "Form1"
   ScaleHeight     =   8595
   ScaleWidth      =   11880
   Begin VB.TextBox Text14 
      Height          =   285
      Left            =   840
      TabIndex        =   2
      Text            =   "10.0"
      Top             =   2040
      Width           =   495
   End
   Begin VB.TextBox Text13 
      Height          =   285
      Left            =   840
      TabIndex        =   38
      Text            =   "0.01"
      Top             =   1680
      Width           =   495
   End
   Begin VB.TextBox Text12 
      Height          =   285
      Left            =   720
      TabIndex        =   36
      Text            =   "0.10"
      Top             =   6720
      Width           =   615
   End
   Begin VB.TextBox Text11 
      Height          =   285
      Left            =   720
      TabIndex        =   1
      Text            =   "5.00"
      Top             =   6360
      Width           =   615
   End
   Begin VB.TextBox Text10 
      Height          =   285
      Left            =   720
      TabIndex        =   34
      Text            =   "0.10"
      Top             =   6120
      Width           =   615
   End
   Begin VB.TextBox Text9 
      Height          =   285
      Left            =   720
      TabIndex        =   33
      Text            =   "0.10"
      Top             =   5880
      Width           =   615
   End
   Begin VB.TextBox Text8 
      Height          =   285
      Left            =   840
      TabIndex        =   29
      Text            =   "0.01"
      Top             =   1320
      Width           =   495
   End
   Begin VB.CommandButton cmdDiffusion 
      Caption         =   "Diffusion"
      Height          =   375
      Left            =   240
      TabIndex        =   27
      Top             =   8160
      Width           =   1095
   End
   Begin VB.TextBox Text7 
      Height          =   285
      Left            =   720
      TabIndex        =   26
      Text            =   "2.00"
      Top             =   7800
      Width           =   615
   End
   Begin VB.TextBox Text6 
      Height          =   285
      Left            =   720
      TabIndex        =   24
      Text            =   "0.20"
      Top             =   7560
      Width           =   615
   End
   Begin VB.TextBox Text5 
      Height          =   285
      Left            =   720
      TabIndex        =   22
      Text            =   "0.01"
      Top             =   7320
      Width           =   615
   End
   Begin VB.TextBox Text4 
      Height          =   285
      Left            =   720
      TabIndex        =   20
      Text            =   "0.10"
      Top             =   7080
      Width           =   615
   End
   Begin VB.TextBox Text2 
      Height          =   285
      Left            =   840
      TabIndex        =   18
      Text            =   "0.50"
      Top             =   600
      Width           =   495
   End
   Begin VB.CommandButton cmdClear 
      Caption         =   "Refresh"
      Height          =   375
      Left            =   240
      TabIndex        =   8
      Top             =   3480
      Width           =   1095
   End
   Begin VB.CommandButton cmdSweep 
      Caption         =   "Sweep"
      Height          =   375
      Left            =   240
      TabIndex        =   9
      Top             =   3120
      Width           =   1095
   End
   Begin VB.CheckBox Check3 
      Caption         =   "Ions"
      Height          =   495
      Left            =   240
      TabIndex        =   16
      Top             =   5280
      Width           =   1335
   End
   Begin VB.CheckBox Check2 
      Caption         =   "Redox"
      Height          =   495
      Left            =   240
      TabIndex        =   15
      Top             =   4800
      Width           =   1095
   End
   Begin VB.CheckBox Check1 
      Caption         =   "Exchange"
      Height          =   495
      Left            =   240
      TabIndex        =   14
      Top             =   4320
      Value           =   1  'Checked
      Width           =   1095
   End
   Begin VB.CommandButton cmdStop 
      Caption         =   "Pause"
      Default         =   -1  'True
      Height          =   375
      Left            =   240
      TabIndex        =   7
      Top             =   3840
      Width           =   1095
   End
   Begin VB.CommandButton cmdStart 
      Caption         =   "Pulse"
      Height          =   375
      Left            =   240
      TabIndex        =   6
      Top             =   2760
      Width           =   1095
   End
   Begin VB.TextBox Text3 
      Height          =   285
      Left            =   840
      TabIndex        =   4
      Text            =   "1.00"
      Top             =   960
      Width           =   495
   End
   Begin VB.TextBox Text1 
      Height          =   285
      Left            =   840
      TabIndex        =   5
      Text            =   "0.45"
      Top             =   240
      Width           =   495
   End
   Begin VB.Label Label7 
      Caption         =   "RFe"
      Height          =   255
      Left            =   240
      TabIndex        =   39
      Top             =   2040
      Width           =   495
   End
   Begin VB.Label Label6 
      Caption         =   "Fe"
      Height          =   255
      Left            =   240
      TabIndex        =   37
      Top             =   1680
      Width           =   615
   End
   Begin VB.Label Label5 
      Caption         =   "oxide"
      Height          =   255
      Left            =   240
      TabIndex        =   35
      Top             =   6720
      Width           =   495
   End
   Begin VB.Label Label4 
      Caption         =   "Rsn"
      Height          =   255
      Left            =   240
      TabIndex        =   32
      Top             =   6360
      Width           =   495
   End
   Begin VB.Label Label3 
      Caption         =   "Rag"
      Height          =   255
      Left            =   240
      TabIndex        =   31
      Top             =   6120
      Width           =   495
   End
   Begin VB.Label Label2 
      Caption         =   "Rag+"
      Height          =   255
      Left            =   240
      TabIndex        =   30
      Top             =   5880
      Width           =   495
   End
   Begin VB.Label Label1 
      Caption         =   "Sigma"
      Height          =   255
      Left            =   240
      TabIndex        =   28
      Top             =   1320
      Width           =   495
   End
   Begin VB.Label lblSnIon 
      Caption         =   "Dsn4"
      Height          =   255
      Left            =   240
      TabIndex        =   25
      Top             =   7800
      Width           =   375
   End
   Begin VB.Label lblGlass 
      Caption         =   "glass"
      Height          =   255
      Left            =   240
      TabIndex        =   23
      Top             =   7560
      Width           =   375
   End
   Begin VB.Label lblAg 
      Caption         =   "Dag"
      Height          =   255
      Left            =   240
      TabIndex        =   21
      Top             =   7320
      Width           =   375
   End
   Begin VB.Label lblAgIon 
      Caption         =   "Dag+"
      Height          =   255
      Left            =   240
      TabIndex        =   19
      Top             =   7080
      Width           =   375
   End
   Begin VB.Label lblAlpha 
      Caption         =   "Alpha"
      BeginProperty Font 
         Name            =   "Times New Roman"
         Size            =   9.75
         Charset         =   255
         Weight          =   400
         Underline       =   0   'False
         Italic          =   0   'False
         Strikethrough   =   0   'False
      EndProperty
      Height          =   255
      Left            =   240
      TabIndex        =   17
      Top             =   600
      Width           =   495
   End
   Begin VB.Label lblX 
      Height          =   495
      Left            =   9600
      TabIndex        =   13
      Top             =   3480
      Width           =   615
   End
   Begin VB.Label lblTime 
      Height          =   495
      Left            =   9600
      TabIndex        =   12
      Top             =   6840
      Width           =   615
   End
   Begin VB.Label lblI 
      Alignment       =   1  'Right Justify
      Height          =   495
      Left            =   1440
      TabIndex        =   11
      Top             =   4320
      Width           =   735
   End
   Begin VB.Label lblC 
      Alignment       =   1  'Right Justify
      Height          =   495
      Left            =   1560
      TabIndex        =   10
      Top             =   240
      Width           =   495
   End
   Begin VB.Label lblLambda 
      Caption         =   "Lambda"
      BeginProperty Font 
         Name            =   "Times New Roman"
         Size            =   9.75
         Charset         =   255
         Weight          =   400
         Underline       =   0   'False
         Italic          =   0   'False
         Strikethrough   =   0   'False
      EndProperty
      Height          =   195
      Left            =   120
      TabIndex        =   3
      Top             =   960
      Width           =   690
   End
   Begin VB.Label lblDo 
      Caption         =   "Dm"
      Height          =   255
      Left            =   240
      TabIndex        =   0
      Top             =   240
      Width           =   255
   End
End
Attribute VB_Name = "Form2"
Attribute VB_GlobalNameSpace = False
Attribute VB_Creatable = False
Attribute VB_PredeclaredId = True
Attribute VB_Exposed = False
    Dim j As Integer
    Dim L As Integer
    Dim m As Integer
    Dim n As Integer
    Dim nmax As Integer
    Dim jmax As Integer
    Dim jint As Integer
    Dim a(1000) As Double
    Dim anew(1000) As Double
    Dim av(1000) As Double
    Dim avnew(1000) As Double
    Dim aa(1000) As Double
    Dim aanew(1000) As Double
    Dim Asum As Double
    Dim Alpha As Double
    Dim b(1000) As Double
    Dim bnew(1000) As Double
    Dim c(1000) As Double
    Dim cnew(1000) As Double
    Dim d(1000) As Double
    Dim dnew(1000) As Double
    Dim Sn(1000) As Double
    Dim Snnew(1000) As Double
    Dim F(1000) As Double
    Dim Fnew(1000) As Double
    Dim Fe(1000) As Double
    Dim Feinit As Double
    Dim Fenew(1000) As Double
    Dim z(1000) As Double
    Dim Chi(1000) As Double
    Dim dChi As Double
    Dim ChiMax As Double
    Dim Tau As Double
    Dim TauBV As Double
    Dim TauStar As Double
    Dim TauFinal As Double
    Dim dTau As Double
    Dim Psi As Double
    Dim PsiExpected As Double
    Dim PsiBV As Double
    Dim PsiV As Double
    Dim X As Double, Y As Double
    Dim La As Double
    Dim Lambda As Double
    Dim Sigma As Double
    Dim Dm As Double
    Dim Da As Double
    Dim DAg As Double
    Dim Dg As Double
    Dim Ds As Double
    Dim Etrans As Double
    Dim Ra As Double
    Dim Ran As Double
    Dim Rc As Double
    Dim RSn As Double
    Dim RFe As Double
    Dim RO As Double
    Dim Rin(1000) As Double
    Dim Rinnew(1000) As Double
    Dim RSi(1000) As Double
    Dim RSinew(1000) As Double
    Dim temp(1000) As Double
    Dim tmp(1000) As Double
    Dim tp(1000) As Double
    Dim Theta As Double
    Dim Pause As Boolean
    Dim Cycle As Boolean
    
    
Private Sub cmdClear_Click()
    

    Scale (-50, 500)-(200, -100)
    
    Line (-5, -100)-(200, 500), &H8000000F, BF
    
    lblC.Caption = ""
    lblX.Caption = ""
    lblI.Caption = ""
    lblTime.Caption = ""
    
    
End Sub


Private Sub cmdStart_Click()
    Dim x1 As Single, y1 As Single
    Dim y2 As Single
    Dim x3 As Single, y3 As Single
    Dim y4 As Single
    Dim y5 As Single
    Dim x6 As Single, y6 As Single
    Dim y7 As Single
    Scale (-50, 500)-(200, -100)
    
    x1 = 0
    y1 = 50
    y2 = 200
    x3 = 150
    y4 = 280
    y5 = 480
    x6 = 150
    y6 = 280
    y7 = -50
    
    DrawWidth = 1

    Line (x1, y7)-(x1, y2)
    Line (x1, y1)-(x3, y1)
    Line (x1, y4)-(x1, y5)
    Line (x1, y4)-(x6, y6)
        
    lblC.Caption = "a(chi)"
    lblX.Caption = "chi"
    lblI.Caption = "psi(tau)"
    lblTime.Caption = "tau"
    
    Lambda = Text3.Text
    If Lambda > 6 Then
        Lambda = 6
    End If
    
    
    DrawWidth = 2
    For j = 2 To jmax Step 2
        X = j
        Y = 280 + 150 * tp(j)
        PSet (X, Y), &H8000000F
    Next
    
    
    DrawWidth = 1
    For j = 1 To jmax - 1
        X = j
        Y = 280 + temp(j) * 150
        PSet (X, Y), &H8000000F
    Next
    
    
    For j = 2 To jmax Step 2
        X = j
        Y = 280 + tmp(j) * 150
        PSet (X, Y), &H8000000F
    Next
    
    
    If Text1.Text < 0.5 Then
        Dm = Text1.Text
    Else
        Dm = 0.49
    End If
    
TauStar = 1#
TauFinal = TauStar * 2#
nmax = 1000
dTau = TauStar / nmax

ChiMax = 6#
dChi = Sqr(dTau / Dm)
jmax = ChiMax / dChi

For j = 0 To jmax
    a(j) = 1#
    Chi(j) = j * dChi
Next
    
    
    
For m = 1 To 2

If Cycle = False Then
    If m = 1 Then
        a(0) = 0#
    Else
        a(0) = 1#
    End If
End If
     
    For n = 1 To nmax
    
        Do
            DoEvents
        Loop While Pause
    
        Tau = n * dTau + (m - 1) * TauStar

If Cycle = True Then
        Theta = Exp(10 - 15 * Tau + 30 * (m - 1) * n * dTau)
        a(0) = Theta / (1 + Theta)
End If

        For j = 1 To jmax - 1
            anew(j) = a(j) + Dm * (a(j + 1) - 2# * a(j) + a(j - 1))
        Next
        
        For j = 1 To jmax - 1
            a(j) = anew(j)
        Next
        
        L = L + 1
        If L = 25 Then
            Psi = (-1.5 * a(0) + 2# * a(1) - 0.5 * a(2)) / dChi
            GetExpected
            
                X = 50 * Tau
                Y = 50 + 15 * PsiExpected
                DrawWidth = 2
                PSet (X, Y), &H8000000F
            
            If Cycle = False Then
                X = 50 * Tau
                Y = 50 + 15 * Psi
            Else
                X = 100 * (Tau + 2 * (1 - m) * n * dTau)
                Y = 50 + 45 * Psi
            End If
            
            DrawWidth = 1
            PSet (X, Y), &H8000000F

            L = 0
        End If
        
    Next
Next
    
    Cycle = False
        
    Cottrell
    
End Sub

Private Sub Cottrell()

    If Text1.Text < 0.5 Then
        Dm = Text1.Text
    Else
        Dm = 0.49
    End If

TauStar = 1#
TauFinal = TauStar * 2#
nmax = 1000
dTau = TauStar / nmax

ChiMax = 6#
dChi = Sqr(dTau / Dm)
jmax = ChiMax / dChi

For j = 0 To jmax
    a(j) = 1#
    Chi(j) = j * dChi
Next

InitPlots
L = 15
Pause = False

For m = 1 To 2

    If m = 1 Then
        a(0) = 0#
    Else
        a(0) = 1#
    End If
    
    For n = 1 To nmax
    
        Do
            DoEvents
        Loop While Pause
    
        Tau = n * dTau + (m - 1) * TauStar
        
        For j = 1 To jmax - 1
            anew(j) = a(j) + Dm * (a(j + 1) - 2# * a(j) + a(j - 1))
        Next
        
        For j = 1 To jmax - 1
            a(j) = anew(j)
        Next
        
        L = L + 1
        If L = 25 Then
            Psi = (-1.5 * a(0) + 2# * a(1) - 0.5 * a(2)) / dChi
            GetExpected
            DoPlots
            L = 0
        End If
        
    Next
Next

End Sub

Public Sub GetExpected()
    pi = 3.14157

    
    If Text1.Text < 0.5 Then
        Dm = Text1.Text
    Else
        Dm = 0.49
    End If
    
    Alpha = Text2.Text
        If Alpha < 0 Then
        Alpha = 0
        End If
        
        If Alpha > 1 Then
        Alpha = 1
        End If
        
    
    If Text3.Text < 6 Then
        If Lambda < 0.00001 Then
        Lambda = 0.00001
        Else
        Lambda = Text3.Text
        End If
    Else
        Lambda = 6
    End If
    
        If Lambda < 0.00001 Then
        Lambda = 0.00001
        End If
    

    If m = 1 Then
        PsiExpected = 1# / Sqr(pi * Tau)
        PsiBV = Lambda * Exp((Lambda ^ 2) * Tau) * erfc(Lambda * Sqr(Tau))
        For j = 0 To jmax Step 2
            z(j) = erf(Chi(j) / (2# * Sqr(Tau)))
            b(j) = 1# - a(j)
        Next
    Else
        PsiExpected = (1# / Sqr(pi * Tau)) - (1# / Sqr(pi * (Tau - TauStar)))
        PsiBV = Lambda * Exp((Lambda ^ 2) * Tau) * erfc(Lambda * Sqr(Tau)) - Lambda * Exp((Lambda ^ 2) * (Tau - TauStar)) * erfc(Lambda * Sqr(Tau - TauStar))

        For j = 0 To jmax Step 2
            z(j) = erf(Chi(j) / (2# * Sqr(Tau))) + erfc(Chi(j) / (2# * Sqr(Tau - TauStar)))
            b(j) = 1# - a(j)
        Next
    End If
               
End Sub

Public Function erf(X As Double) As Double

    Dim z As Double
    
    z = 1# / (1# + 0.47047 * X)
    erf = 0.3480242 * z - 0.0958798 * (z ^ 2#) + 0.7478556 * (z ^ 3#)
    erf = 1# - erf * Exp(-(X ^ 2#))
    
End Function

Public Function erfc(X As Double) As Double

    erfc = 1# - erf(X)
    
End Function

Private Sub InitPlots()
    
        For j = 1 To jmax - 1
            temp(j) = 1
        Next
        
        For j = 1 To jmax - 1
            tmp(j) = 0
        Next
        
        For j = 2 To jmax Step 2
            tp(j) = 1
        Next
        
        
        
End Sub

Private Sub DoPlots()
    
    If Check2 = 1 Then
    
    DrawWidth = 2
    For j = 2 To jmax Step 2
        X = j
        Y = 280 + 150 * tp(j)
        PSet (X, Y), &H8000000F
    Next

            For j = 0 To jmax Step 2
                tp(j) = z(j)
            Next

        For j = 2 To jmax Step 2
            X = j
            Y = 280 + 150 * z(j)
            PSet (X, Y), vbYellow
        Next
        
    End If

    DrawWidth = 1
    For j = 1 To jmax - 1
        X = j
        Y = 280 + temp(j) * 150
        PSet (X, Y), &H8000000F
    Next

            For j = 1 To jmax - 1
                temp(j) = a(j)
            Next

    If Check1 = 1 Then
        For j = 1 To jmax - 1
            X = j
            Y = 280 + a(j) * 150
            PSet (X, Y)
        Next
    End If

    For j = 2 To jmax Step 2
        X = j
        Y = 280 + tmp(j) * 150
        PSet (X, Y), &H8000000F
    Next

    
    If Check1 = 1 Then
        For j = 0 To jmax Step 2
            X = j
            Y = 280
            PSet (X, Y)
        Next
    End If
    
            For j = 2 To jmax Step 2
                tmp(j) = b(j)
            Next
            
    If Check1 = 1 Then
        For j = 2 To jmax Step 2
            X = j
            Y = 280 + b(j) * 150
            PSet (X, Y), vbBlue
        Next
    End If



        If Check2 = 1 Then
            X = 50 * Tau
            Y = 50 + 15 * PsiExpected
            DrawWidth = 2
            PSet (X, Y), vbYellow
        End If
        
        
        If Check3 = 1 Then
            If m = 1 Then
                X = 50 * Tau
                Y = 50 + 15 * PsiBV
                DrawWidth = 2
                PSet (X, Y), vbRed
            End If
        End If
        

        If Check1 = 1 Then
            DrawWidth = 1
            
            If Cycle = False Then
                X = 50 * Tau
                Y = 50 + 15 * Psi
            Else
                X = 100 * (Tau + 2 * (1 - m) * n * dTau)
                Y = 50 + 45 * Psi
            End If
            
            PSet (X, Y)
        End If
        
    
End Sub

Private Sub cmdStop_Click()
    Pause = Not Pause
End Sub

Private Sub BV()
    E0 = 10#
    Slope = 15#
    
    If Text1.Text < 0.5 Then
        Dm = Text1.Text
    Else
        Dm = 0.49
    End If
    
    Alpha = Text2.Text
        If Alpha < 0 Then
        Alpha = 0
        End If
        
        If Alpha > 1 Then
        Alpha = 1
        End If
        
    Lambda = Text3.Text
        If Lambda < 0.00001 Then
        Lambda = 0.00001
        End If
    

    
TauStar = 1#
TauFinal = TauStar * 2#
nmax = 1000
dTau = TauStar / nmax

ChiMax = 6#
dChi = Sqr(dTau / Dm)
jmax = ChiMax / dChi

For j = 0 To jmax
    a(j) = 1#
    av(j) = 1#
    Chi(j) = j * dChi
Next

InitPlots
L = 15
Pause = False

For m = 1 To 2
   
    For n = 1 To nmax
    
        Do
            DoEvents
        Loop While Pause
    
        TauBV = n * dTau + (m - 1) * TauStar
        Theta = Exp(E0 - Slope * TauBV + 2 * Slope * (m - 1) * n * dTau)
        La = Lambda * (1 + Theta) / (Theta ^ Alpha)
        Psi = (-1.5 * (Theta / (1 + Theta)) + 2# * a(1) - 0.5 * a(2)) / (dChi * (1 + 1.5 / (dChi * La)))
        
        av(0) = Theta / (1 + Theta)
        
        For j = 1 To jmax - 1
            avnew(j) = av(j) + Dm * (av(j + 1) - 2# * av(j) + av(j - 1))
        Next
        
        For j = 1 To jmax - 1
            av(j) = avnew(j)
        Next
        
            a(0) = (Psi / La) + Theta / (1 + Theta)
        For j = 1 To jmax - 1
            anew(j) = a(j) + Dm * (a(j + 1) - 2# * a(j) + a(j - 1))
        Next
        
        For j = 1 To jmax - 1
            a(j) = anew(j)
        Next
        
        For j = 0 To jmax
        
            If a(j) < 0 Then
            a(j) = 0
            End If
            If a(j) > 1 Then
            a(j) = 1
            End If

        Next

        
        L = L + 1
        If L = 25 Then
            PsiV = (-1.5 * av(0) + 2# * av(1) - 0.5 * av(2)) / dChi
            
            Plots
            L = 0
        End If
        
    Next
Next

End Sub

Private Sub Plots()
    
    If Check1 = 1 Then
    
    DrawWidth = 1

    For j = 1 To jmax - 1
        X = j
        Y = 280 + temp(j) * 150
        PSet (X, Y), &H8000000F
    Next

    For j = 1 To jmax - 1
        temp(j) = av(j)
    Next

    For j = 1 To jmax - 1
        X = j
        Y = 280 + av(j) * 150
        PSet (X, Y)
    Next
        
            X = 100 * (TauBV + 2 * (1 - m) * n * dTau)
            Y = 50 + 45 * PsiV
            PSet (X, Y)


    End If

    
    If Check3 = 1 Then
    
    DrawWidth = 1

    For j = 2 To jmax - 2
        X = j
        Y = 280 + tp(j) * 150
        PSet (X, Y), &H8000000F
    Next

            For j = 2 To jmax - 2
                tp(j) = a(j)
            Next

        For j = 2 To jmax - 2 Step 2
            X = j
            Y = 280 + a(j) * 150
            PSet (X, Y), vbRed
        Next
        
            X = 100 * (TauBV + 2 * (1 - m) * n * dTau)
            Y = 50 + 45 * Psi
            PSet (X, Y), vbRed

    End If

        
    
End Sub

Private Sub cmdSweep_Click()

    Dim x1 As Single, y1 As Single
    Dim y2 As Single
    Dim x3 As Single, y3 As Single
    Dim y4 As Single
    Dim y5 As Single
    Dim x6 As Single, y6 As Single
    Dim y7 As Single
    Scale (-50, 500)-(200, -100)
    
    x1 = 0
    y1 = 50
    y2 = 200
    x3 = 150
    y4 = 280
    y5 = 480
    x6 = 150
    y6 = 280
    y7 = -50
    
    DrawWidth = 1
  
    Line (x1, y7)-(x1, y2)
    Line (x1, y1)-(x3, y1)
    Line (x1, y4)-(x1, y5)
    Line (x1, y4)-(x6, y6)
        
    lblC.Caption = "a(chi)"
    lblX.Caption = "chi"
    lblI.Caption = "psi(theta)"
    lblTime.Caption = "theta"
    
    Check2 = 0
    
    
    DrawWidth = 2
    For j = 2 To jmax Step 2
        X = j
        Y = 280 + 150 * tp(j)
        PSet (X, Y), &H8000000F
    Next
    
    
    DrawWidth = 1
    For j = 1 To jmax - 1
        X = j
        Y = 280 + temp(j) * 150
        PSet (X, Y), &H8000000F
    Next
    
    
    For j = 1 To jmax - 1
        X = j
        Y = 280 + tmp(j) * 150
        PSet (X, Y), &H8000000F
    Next
    
    
    If Text1.Text < 0.5 Then
        Dm = Text1.Text
    Else
        Dm = 0.49
    End If
    
    Lambda = Text3.Text
    If Lambda > 6 Then
        Lambda = 6
    End If
    

    
TauStar = 1#
TauFinal = TauStar * 2#
nmax = 1000
dTau = TauStar / nmax

ChiMax = 6#
dChi = Sqr(dTau / Dm)
jmax = ChiMax / dChi

For j = 0 To jmax
    a(j) = 1#
    Chi(j) = j * dChi
Next
    
L = 15

For m = 1 To 2

If Cycle = False Then
    If m = 1 Then
        a(0) = 0#
    Else
        a(0) = 1#
    End If
End If
    
    For n = 1 To nmax
    
        Do
            DoEvents
        Loop While Pause
    
        Tau = n * dTau + (m - 1) * TauStar

If Cycle = True Then
        Theta = Exp(10 - 15 * Tau + 30 * (m - 1) * n * dTau)
        a(0) = Theta / (1 + Theta)
End If

        For j = 1 To jmax - 1
            anew(j) = a(j) + Dm * (a(j + 1) - 2# * a(j) + a(j - 1))
        Next
        
        For j = 1 To jmax - 1
            a(j) = anew(j)
        Next
        
        L = L + 1
        If L = 25 Then
            Psi = (-1.5 * a(0) + 2# * a(1) - 0.5 * a(2)) / dChi
            GetExpected
                X = 50 * Tau
                Y = 50 + 15 * PsiExpected
                DrawWidth = 2
                PSet (X, Y), &H8000000F

                
            If Cycle = False Then
                X = 50 * Tau
                Y = 50 + 15 * Psi
            Else
                X = 100 * (Tau + 2 * (1 - m) * n * dTau)
                Y = 50 + 45 * Psi
            End If
            
            DrawWidth = 1
            PSet (X, Y), &H8000000F

            L = 0
        End If
        
    Next
Next
    
    Cycle = True

    BV
    

End Sub

Private Sub cmdDiffusion_Click()
    Dim x1 As Single, y1 As Single
    Dim y2 As Single
    Dim x3 As Single, y3 As Single
    Dim y4 As Single
    Dim y5 As Single
    Dim x6 As Single, y6 As Single
    Dim y7 As Single
    Dim x8 As Single, y8 As Single
    Scale (-50, 500)-(200, -100)
    
    x1 = 0
    y1 = 50
    y2 = 200
    x3 = 150
    y4 = 280
    y5 = 480
    x6 = 150
    y6 = 280
    y7 = -50
    
    DrawWidth = 1

    Line (x1, y7)-(x1, y2)
    Line (x1, y1)-(x3, y1)
    Line (x1, y4)-(x1, y5)
    Line (x1, y4)-(x6, y6)
            
    lblC.Caption = "a(chi)"
    lblX.Caption = "chi"
    lblI.Caption = "psi(tau)"
    lblTime.Caption = "tau"
    
    
    
    Lambda = Text3.Text
    If Lambda > 6 Then
        Lambda = 6
    End If
     
    If Text1.Text < 0.5 Then
        Dm = Text1.Text
    Else
        Dm = 0.49
    End If
       
  
    
    Cycle = False
        
        
    Diffusion
    
End Sub

Private Sub Diffusion()

    If Text1.Text < 0.5 Then
        Dm = Text1.Text
    Else
        Dm = 0.49
    End If
    
    If Text4.Text < 1# Then
        Da = Text4.Text
    Else
        Da = 0.99
    End If
    
    If Text5.Text < 1# Then
        DAg = Text5.Text
    Else
        DAg = 0.99
    End If
    
    If Text6.Text < 1# Then
        Dg = Text6.Text
    Else
        Dg = 0.99
    End If
    
    If Text7.Text < 10# Then
        Ds = Text7.Text
    Else
        Ds = 9.9
    End If
    
    If Text8.Text < 10# Then
        Sigma = Text8.Text
    Else
        Sigma = 9.9
    End If
    
    If Text9.Text < 1# Then
        Ra = Text9.Text * 10000
    Else
        Ra = 0.99 * 10000
    End If
    
    If Text10.Text < 1# Then
        Ran = Text10.Text * 10000
    Else
        Ran = 0.99 * 10000
    End If
    
    If Text11.Text < 100# Then
        RSn = Text11.Text * 10000
    Else
        RSn = 99 * 10000
    End If
    
    If Text12.Text < 0.5 Then
        cmax = Text12.Text
    Else
        cmax = 0.49
    End If

    
TauStar = 1#
nmax = 4000
dTau = TauStar / nmax

ChiMax = 6#
ChiStar = 1#
dChi = Sqr(dTau / Dm)
jmax = ChiMax / dChi
jint = ChiStar / dChi
dmax = 0.16

    
    If Text13.Text < 0.1 Then
        Feinit = Text13.Text
    Else
        Feinit = 0.09
    End If
    
    If Text14.Text < 100# Then
        RFe = Text14.Text * 10000
    Else
        RFe = 99 * 10000
    End If

    x8 = jint / 2#
    y4 = 280
    y8 = 380
    
    Line (x8, y4)-(x8, y8)

For j = 0 To jmax
    a(j) = 0#
    aa(j) = 0#
    c(j) = 0#
    Chi(j) = j * dChi
Next

aa(0) = 1#

For j = 0 To jint
    b(j) = 0#
    d(j) = 0#
    Sn(j) = 0#
    F(j) = 0#
    Fe(j) = 0#
    Rin(j) = 0#
    RSi(j) = 0#
Next

For j = jint + 1 To jmax
    Sn(j) = 0.01
    b(j) = dmax - 4 * Sn(j)
    d(j) = dmax
    F(j) = 0.001
    Fe(j) = Feinit
    Rin(j) = 0#
    RSi(j) = 0#
Next
    
            b(jmax) = dmax
            b(jmax - 1) = dmax - 2 * Sn(j)
      

InitPlots
m = 1
L = 15
Pause = False

    For n = 1 To nmax
    
        Do
            DoEvents
        Loop While Pause
    
        Tau = n * dTau
        Asum = 0
        
        For j = 1 To jint - 1
            Rinnew(j) = Ra * aa(j) * (cmax - c(j)) - Ran * a(j) * c(j)
            Asum = Asum + a(j) + b(j) - c(j) + (4 * Sn(j))
            aanew(j) = aa(j) + DAg * Dm * (aa(j + 1) - 2# * aa(j) + aa(j - 1)) - dTau * Rin(j)
            anew(j) = a(j) + Da * Dm * ((a(j + 1) - 2# * a(j) + a(j - 1)) - Sigma * (a(j) * (a(j) + b(j) - c(j) + (4 * Sn(j))) - Asum * (a(j + 1) - a(j - 1)) / 2)) + dTau * Rin(j) + 4 * dTau * RSn * aa(j) * Sn(j)
            bnew(j) = b(j) + Dm * ((b(j + 1) - 2# * b(j) + b(j - 1)) - Sigma * (b(j) * (a(j) + b(j) - c(j) + (4 * Sn(j))) - Asum * (b(j + 1) - b(j - 1)) / 2))
            Snnew(j) = Sn(j) + DSn * Dm * ((Sn(j + 1) - 2# * Sn(j) + Sn(j - 1)) - 4 * Sigma * (Sn(j) * (a(j) + b(j) - c(j) + (4 * Sn(j))) - Asum * (Sn(j + 1) - Sn(j - 1)) / 2)) - dTau * RSn * aa(j) * Sn(j)
            cnew(j) = c(j) + dTau * Rin(j)
        Next
        
        For j = jint To jmax - 1
            RSinew(j) = RFe * Fe(j) * a(j) * (d(j) - b(j) - 4 * Sn(j))
            Asum = Asum + a(j) + b(j) - d(j) + (4 * Sn(j))
            aanew(j) = aa(j) + DAg * Dm * (aa(j + 1) - 2# * aa(j) + aa(j - 1)) + dTau * RSi(j)
            anew(j) = a(j) + Da * Dm * ((a(j + 1) - 2# * a(j) + a(j - 1)) - Sigma * (a(j) * (a(j) + b(j) - d(j) + (4 * Sn(j))) - Asum * (a(j + 1) - a(j - 1)) / 2)) + 4 * dTau * RSn * aa(j) * Sn(j) - dTau * RSi(j)
            bnew(j) = b(j) + Dm * ((b(j + 1) - 2# * b(j) + b(j - 1)) - Sigma * (b(j) * (a(j) + b(j) - d(j) + (4 * Sn(j))) - Asum * (b(j + 1) - b(j - 1)) / 2))
            Snnew(j) = Sn(j) + DSn * Dm * ((Sn(j + 1) - 2# * Sn(j) + Sn(j - 1)) - 4 * Sigma * (Sn(j) * (a(j) + b(j) - d(j) + (4 * Sn(j))) - Asum * (Sn(j + 1) - Sn(j - 1)) / 2)) - dTau * RSn * aa(j) * Sn(j)
            dnew(j) = d(j) - dTau * RSi(j)
        Next
        
        For j = 1 To jmax - 1
        If anew(j) < 0 Then
            a(j) = 0
        Else
            a(j) = anew(j)
        End If
        If aanew(j) < 0 Then
            aa(j) = 0
        Else
            aa(j) = aanew(j)
        End If
        If bnew(j) < 0 Then
            b(j) = 0
        Else
            b(j) = bnew(j)
        End If
            b(jmax) = dmax
            b(jmax - 1) = dmax - 2 * Sn(j)
        If cnew(j) < 0 Then
            c(j) = 0
        Else
            If cnew(j) > cmax Then
            c(j) = cmax
            Else
            c(j) = cnew(j)
            End If
        End If
        If dnew(j) < 0 Then
            d(j) = 0
        Else
            If dnew(j) > dmax Then
            d(j) = dmax
            Else
            d(j) = dnew(j)
            End If
        End If
        If Snnew(j) < 0 Then
            Sn(j) = 0
        Else
            Sn(j) = Snnew(j)
        End If
        If RSinew(j) < 0 Then
            RSi(j) = 0
        Else
            RSi(j) = RSinew(j)
        End If
        
        Next
        
        L = L + 1
        
        If L = 25 Then
            Etrans = 0
            Psi = 0
        
            For j = 1 To jint - 1
                Etrans = Etrans + dTau * Rin(j) + 4 * dTau * RSn * aa(j) * Sn(j)
            Next
            
            For j = jint To jmax - 1
                Psi = Psi + dTau * RFe * Fe(j) * a(j) * (d(j) - b(j) - 4 * Sn(j))
            Next
        
            DoPlot
            L = 0
            
        End If
        
    Next


End Sub


Private Sub DoPlot()
    
    
    DrawWidth = 1
    
    Scale (-50, 500)-(200, -100)
    
    x1 = 0
    y4 = 280
    x6 = 150
    y6 = 280
    
    DrawWidth = 1
    Line (x1, y4)-(x6, y6)
  
    DrawWidth = 2
                X = 50 * Tau
                Y = 50 + 15000 * Psi
                PSet (X, Y), vbBlue

                X = 50 * Tau
                Y = 50 + 7500 * Etrans
                PSet (X, Y), vbRed
        
    Scale (-50, 500)-(200, -100)
    
    x1 = 0
    y4 = 280
    x6 = 150
    y6 = 280
    
    DrawWidth = 1
    Line (x1, y4)-(x6, y6)
  
        
        
If Check1 = 1 Then
    
    DrawWidth = 2
    
    For j = 1 To jmax - 1
        X = j / 2#
        Y = 280 + temp(j) * 150
        PSet (X, Y), &H8000000F
    Next
    
    For j = 2 To jmax - 1
        X = (j - 0.5) / 2#
        Y = 280 + 75 * (temp(j - 1) + temp(j))
        PSet (X, Y), &H8000000F
    Next
    
    For j = 1 To jmax - 1
        X = j / 2#
        Y = 280 + aa(j) * 150
        PSet (X, Y), vbBlue
    Next

    For j = 2 To jmax - 1
        X = (j - 0.5) / 2#
        Y = 280 + 75 * (aa(j - 1) + aa(j))
        PSet (X, Y), vbBlue
    Next

            For j = 1 To jmax - 1
                temp(j) = aa(j)
            Next

    DrawWidth = 2
    
    For j = 1 To jmax - 1
        X = j / 2#
        Y = 280 + tmp(j) * 150
        PSet (X, Y), &H8000000F
    Next
    
    For j = 2 To jmax - 1
        X = (j - 0.5) / 2#
        Y = 280 + 75 * (tmp(j - 1) + tmp(j))
        PSet (X, Y), &H8000000F
    Next
    
    For j = 1 To jmax - 1
        X = j / 2#
        Y = 280 + b(j) * 150
        PSet (X, Y), vbGreen
    Next

    For j = 2 To jmax - 1
        X = (j - 0.5) / 2#
        Y = 280 + 75 * (b(j - 1) + b(j))
        PSet (X, Y), vbGreen
    Next

            For j = 1 To jmax - 1
                tmp(j) = b(j)
            Next

    DrawWidth = 2
    
    For j = 1 To jmax - 1
        X = j / 2#
        Y = 280 + tp(j) * 150
        PSet (X, Y), &H8000000F
    Next
    
    For j = 2 To jmax - 1
        X = (j - 0.5) / 2#
        Y = 280 + 75 * (tp(j - 1) + tp(j))
        PSet (X, Y), &H8000000F
    Next
    
  
    For j = 1 To jmax - 1
        X = j / 2#
        Y = 280 + a(j) * 150
        PSet (X, Y), vbRed
    Next

    For j = 2 To jmax - 1
        X = (j - 0.5) / 2#
        Y = 280 + 75 * (a(j - 1) + a(j))
        PSet (X, Y), vbRed
    Next

            For j = 1 To jmax - 1
                tp(j) = a(j)
            Next

Else
End If
    
    
    
If Check2 = 1 Then
    
    DrawWidth = 2
    
    For j = 1 To jmax - 1
        X = j / 2#
        Y = 280 + 150 * temp(j) / 100
        PSet (X, Y), &H8000000F
    Next

    For j = 1 To jmax - 1
        X = (j - 0.5) / 2#
        Y = 280 + 75 * (temp(j - 1) + temp(j)) / 100
        PSet (X, Y), &H8000000F
    Next
        
        For j = 1 To jint - 1
            X = j / 2#
            Y = 280 + 150 * Ra * cmax / 100
            PSet (X, Y), vbGreen
        Next
        
        For j = 2 To jint - 1
            X = (j - 0.5) / 2#
            Y = 280 + 75 * (Ra * cmax + Ra * cmax) / 100
            PSet (X, Y), vbGreen
        Next
        
        For j = jint To jmax - 1
            X = j / 2#
            Y = 280 + 150 * (4 * RSn * aa(j) * Sn(j) - 0 * RSi(j)) * 1000
            PSet (X, Y), vbGreen
        Next
        
        For j = 1 To jint - 1
            temp(j) = Ra * cmax
        Next

        For j = jint To jmax - 1
            temp(j) = (4 * RSn * aa(j) * Sn(j) - 0 * RSi(j)) * 100000
        Next


    DrawWidth = 2
    
    For j = 1 To jmax - 1
        X = j / 2#
        Y = 280 + 60 * tmp(j)
        PSet (X, Y), &H8000000F
    Next

    For j = 2 To jmax - 1
        X = (j - 0.5) / 2#
        Y = 280 + 30 * (tmp(j - 1) + tmp(j))
        PSet (X, Y), &H8000000F
    Next

        For j = 1 To jmax - 1
            X = j / 2#
            Y = 280 + 60 * RSn * aa(j) * Sn(j) + 100000 * Rin(j)
            PSet (X, Y), vbRed
        Next
        
        For j = 2 To jmax - 1
            X = (j - 0.5) / 2#
            Y = 280 + 30 * RSn * ((aa(j - 1) * Sn(j - 1)) + (aa(j) * Sn(j))) + 50000 * Rin(j - 1) + 50000 * Rin(j)
            PSet (X, Y), vbRed
        Next
        
        
        For j = 1 To jmax - 1
            tmp(j) = RSn * aa(j) * Sn(j) + (100000 * Rin(j) / 60)
        Next

Else
End If
    
    
    
If Check3 = 1 Then
    
    DrawWidth = 2
    
    For j = 2 To jmax - 1 Step 2
        X = j / 2#
        Y = 280 + 750 * temp(j)
        PSet (X, Y), &H8000000F
    Next

        For j = 2 To jmax - 1 Step 2
            X = j / 2#
            Y = 280 + 750 * a(j)
            PSet (X, Y), vbRed
        Next
        
        For j = 2 To jmax - 1 Step 2
            temp(j) = a(j)
        Next
    
    DrawWidth = 2
    
    For j = 2 To jmax - 1 Step 2
        X = j / 2#
        Y = 280 + 750 * tmp(j)
        PSet (X, Y), &H8000000F
    Next

        For j = 2 To jmax - 1 Step 2
            X = j / 2#
            Y = 280 + 750 * b(j)
            PSet (X, Y), vbBlue
        Next
        
        For j = 2 To jmax - 1 Step 2
            tmp(j) = b(j)
        Next

    DrawWidth = 2
    For j = 1 To jmax - 1
        X = j / 2#
        Y = 280 + 1500 * tp(j)
        PSet (X, Y), &H8000000F
    Next

    For j = 2 To jmax - 1
        X = (j - 0.5) / 2#
        Y = 280 + 750 * (tp(j - 1) + tp(j))
        PSet (X, Y), &H8000000F
    Next

        For j = 1 To jmax - 1
            X = j / 2#
            Y = 280 + 1500 * Sn(j)
            PSet (X, Y), vbYellow
        Next
        
        For j = 2 To jmax - 1
            X = (j - 0.5) / 2#
            Y = 280 + 750 * (Sn(j - 1) + Sn(j))
            PSet (X, Y), vbYellow
        Next
        
        For j = 1 To jmax - 1
            tp(j) = Sn(j)
        Next

Else
End If
    
    
End Sub