Public Structure State
        Dim posX As Single
        Dim posY As Single
        Dim velX As Single
        Dim velY As Single
    End Structure

    Private Structure Derivative
        Dim dx As Single ' dx/dt = velocity
        Dim dy As Single
        Dim dvx As Single ' dv/dt = acceleration
        Dim dvy As Single
    End Structure

    Private Function evaluate(ByVal initial As State, ByVal t As Single, ByVal dt As Single, ByVal d As Derivative) As Derivative
        Dim state As New State()
        state.posX = initial.posX + d.dx * dt
        state.posY = initial.posY + d.dy * dt
        state.velX = initial.velX + d.dvx * dt
        state.velY = initial.velY + d.dvy * dt

        Dim output As New Derivative()
        output.dx = state.velX
        output.dy = state.velY
        output.dvy = accelerationY(state, t + dt)
        output.dvx = accelerationX(state, t + dt)
        Return output
    End Function

    Private Function accelerationX(ByVal state As State, ByVal t As Single) As Single
        Dim k As Single = 0 'Ball.Force.Net.X '
        Const b As Single = 1
        Return (k * state.posX) - b * (state.velX)
    End Function

    Private Function accelerationY(ByVal state As State, ByVal t As Single) As Single
        Dim k As Single = Calculate_Gravitational_Force(Ball.Mass, EARTH_GRAVITY) '- Calculate_Air_Resistance(Ball.Drag_Coefficient, AIR_DENSITY, Ball.Area, Ball.Velocity.Y) '10
        Const b As Single = 1
        Return (k * state.posY) - b * (state.velY)
    End Function

    Private Function integrate(ByVal state As State, ByVal t As Single, ByVal dt As Single) As State
        Dim a As Derivative
        Dim b As Derivative
        Dim c As Derivative
        Dim d As Derivative

        a = (evaluate(state, t, 0.0F, New Derivative()))
        b = (evaluate(state, t, dt * 0.5F, a))
        c = (evaluate(state, t, dt * 0.5F, b))
        d = (evaluate(state, t, dt, c))

        Dim dxdt As Single = (1.0F / 6.0F) * (a.dx + 2.0F * (b.dx + c.dx) + d.dx)
        Dim dydt As Single = (1.0F / 6.0F) * (a.dy + 2.0F * (b.dy + c.dy) + d.dy)

        Dim dvxdt As Single = (1.0F / 6.0F) * (a.dvx + 2.0F * (b.dvx + c.dvx) + d.dvx)
        Dim dvydt As Single = (1.0F / 6.0F) * (a.dvy + 2.0F * (b.dvy + c.dvy) + d.dvy)

        state.posX += dxdt * dt
        state.posY += dydt * dt
        state.velX += dvxdt * dt
        state.velY += dvydt * dt

        Return state

    End Function