Option Explicit

Dim q(129) As Double, p(129) As Double
Dim a(129) As Double, ad(129) As Double
Dim E(129) As Double
Dim n As Integer

Private Sub CommandButton1_Click()
    Dim w(7) As Double
    Dim c(16) As Double
    Dim d(16) As Double
    Dim dt As Double
    Dim i As Long, k As Integer, j As Integer
    Dim pi As Double
    Dim ang As Double
    Dim ii As Integer
    
    Application.ScreenUpdating = False
    
    pi = 3.14159265358979
    
    w(1) = 0.311790812418427
    w(2) = -1.55946803821447
    w(3) = -1.6789692825964
    w(4) = 1.66335809963315
    w(5) = -1.06458714789183
    w(6) = 1.36934946416871
    w(7) = 0.629030650210433
    w(0) = 1# - 2# * (w(1) + w(2) + w(3) + w(4) + w(5) + w(6) + w(7))
    
    d(1) = w(7)
    d(15) = w(7)
    d(2) = w(6)
    d(14) = w(6)
    d(3) = w(5)
    d(13) = w(5)
    d(4) = w(4)
    d(12) = w(4)
    d(5) = w(3)
    d(11) = w(3)
    d(6) = w(2)
    d(10) = w(2)
    d(7) = w(1)
    d(9) = w(1)
    d(8) = w(0)
    d(16) = 0#
    
    c(1) = 0.5 * w(7)
    c(16) = 0.5 * w(7)
    c(2) = 0.5 * (w(7) + w(6))
    c(15) = 0.5 * (w(7) + w(6))
    c(3) = 0.5 * (w(6) + w(5))
    c(14) = 0.5 * (w(6) + w(5))
    c(4) = 0.5 * (w(5) + w(4))
    c(13) = 0.5 * (w(5) + w(4))
    c(5) = 0.5 * (w(4) + w(3))
    c(12) = 0.5 * (w(4) + w(3))
    c(6) = 0.5 * (w(3) + w(2))
    c(11) = 0.5 * (w(3) + w(2))
    c(7) = 0.5 * (w(2) + w(1))
    c(10) = 0.5 * (w(2) + w(1))
    c(8) = 0.5 * (w(1) + w(0))
    c(9) = 0.5 * (w(1) + w(0))
    
    n = 32
    
    q(0) = 0
    q(n + 1) = 0
    
    For j = 1 To n
        q(j) = Sin(pi * j / (n + 1))
        p(j) = 0#
    Next j
    
    dt = Sqr(1 / 8)
    ii = 1
    
    For i = 0 To 30000

        For k = 1 To 16
            For j = 1 To n
                q(j) = q(j) + dt * c(k) * p(j)
            Next j
            For j = 1 To n
                p(j) = p(j) + dt * d(k) * f(j)
            Next j
        Next k

        
        For k = 1 To n
            a(k) = 0
            ad(k) = 0
        Next k
        
        For k = 1 To n
            For j = 1 To n
                a(k) = a(k) + q(j) * Sin(j * k * pi / (n + 1))
                ad(k) = ad(k) + p(j) * Sin(j * k * pi / (n + 1))
            Next j
            E(k) = Sqr(2 / (n + 1)) * ((ad(k) ^ 2) / 2 + 2 * (a(k) ^ 2) * (Sin(pi * k / 2 / (n + 1)) ^ 2))
        Next k
        
        If i Mod 10 = 0 Then
            Worksheets("Sheet1").Cells(ii + 1, 2) = dt * i
                For k = 1 To 10
                Worksheets("Sheet1").Cells(ii + 1, 2 + k) = E(k)
            Next k
        ii = ii + 1
        End If

    Next i


    Application.ScreenUpdating = True
    
End Sub

Function f(k As Integer) As Double
    f = q(k + 1) - 2# * q(k) + q(k - 1) + 0.25 * ((q(k + 1) - q(k)) ^ 2 - (q(k) - q(k - 1)) ^ 2)
End Function