Function RenewalFunction(alpha As Single, lambda As Single, t As Single) As Single
Dim dt As Single
Dim i As Integer
Dim mttf As Single
Const MaxDim As Integer = 200
Const eps As Single = 0.001
Dim cnt As Integer
Dim W(0 To MaxDim) As Single
Dim f(0 To MaxDim) As Single
Dim F1(0 To MaxDim) As Single
Dim prev As Single
Dim w0 As Single
cnt = 0
dt = t / MaxDim
If t <= 0# Then
    RenewalFunction = 0#
    Exit Function
End If

mttf = Gamma(1# / alpha + 1#) / lambda
' Calculate initial W()
For i = 0 To MaxDim
    W(i) = (Gamma(1# / alpha + 1#) * dt * i / mttf) ^ alpha
    F1(i) = 1# - Exp(-(lambda * i * dt) ^ alpha)
    f(i) = alpha * lambda * (lambda * i * dt) ^ (alpha - 1) * (1 - F1(i))
Next
prev = W(MaxDim)
Do
    cnt = cnt + 1
    prev = W(MaxDim)
    For i = MaxDim To 1 Step -1
        W(i) = F1(i) + Conv(i, 0, i, dt, W, f)
    Next
Loop While Abs(prev - W(MaxDim)) / prev > eps Or cnt < 2
RenewalFunction = W(MaxDim)
End Function

Function Conv(t As Integer, a As Integer, b As Integer, h As Single, f() As Single, g() As Single)
' Calculate \int_a^b f(t-x)g(x)dx
' f() and g() are arrays
' a = integer value for pointer for lower limit, i.e., g(a)
' b = integer value for pointer for upper limit, i.e., g(b)
' t = integer value for pointer to t
' h = step-length in f() and g(), i.e., (real) distance between points
Dim i As Integer
Dim j As Integer
Dim sM As Single
Dim a1 As Single
Dim a2 As Single
Dim b1 As Single
Dim b2 As Single
Dim hinv As Single
hinv = 1# / h
sM = 0#
For j = a + 1 To b
    b1 = f(t - j + 1)
    a1 = (f(t - j) - f(t - j + 1)) * hinv
    b2 = g(j - 1)
    a2 = (g(j) - g(j - 1)) * hinv
    sM = sM + h * (h * (0.3333333 * a1 * a2 * h + 0.5 * (a1 * b2 + b1 * a2)) + b1 * b2)
Next
Conv = sM
End Function