Practical implementation of Vellekoop-Nieuwenhuis model/interpolation

Question

Have read the 2006 VELLEKOOP-NIEUWENHUIS paper (Efficient Pricing of Derivatives on Assets with Discrete Dividends) (Download) many times re Discrete dividends on American Options, but remain baffled as to how to implement this with code. The verbose paragraph in the paper about implementation is in image below.

Read most of the citations to look for some hints to no avail, but one source did provide a formula for one possible interpolation without saying what k,j, V or any of the variables in it actually are. See Image below.

One blogger tried to do it, but it is flawed and as cryptic as the paper.

Have included code for CRR Binomial Tree for American Options below, looking to implement that cash div methodology at the appropriate step indicated in accordance with the VELLEKOOP interpolation scheme.

Please let me know the correct way to do this if you know, because I am stuck. Thanks.

import numpy as np

AmeEurFlag='a'
CallPutFlag='P' 
S = 100 
X = 100 
T = 1.0
r = 0.05 
c = r
v = 0.2
n = 10
dt = T / n
#Assume one discrete cash div of $2 at Step 5 (counting from 0)
Div=2
DivStep= 4*dt   


n_list = np.arange(0, (n + 1), 1)
if CallPutFlag == 'C':
    z = 1
elif CallPutFlag == 'P':
    z = -1

# The up and down factors for S 
u = np.exp(v*np.sqrt(dt))
d = 1./u
p = (np.exp((c)*dt)-d) / (u-d) 
df = np.exp(-r * dt)





# Final Columns Of Trees
max_pay_off_list = []

for i in n_list:
    i = i.astype('int')
    max_pay_off = np.maximum(0, z * (S * u ** i * d ** (n - i) - X))
    max_pay_off_list.append(max_pay_off)



# Backwards Recrusion
for j in np.arange(n - 1, 0 - 1, -1):
    
    if j==4: print('\n',"Some Interpolation Magic should have happened here at Step ", i) 
    for i in np.arange(0, j + 1, 1):
        i = i.astype(int)  # Need to be converted to a integer
        if AmeEurFlag == 'e':
            max_pay_off_list[i] = (p * max_pay_off_list[i + 1] + (1 - p) * max_pay_off_list[i]) * df
        elif AmeEurFlag == 'a':
            max_pay_off_list[i] = np.maximum((z * (S * u ** i * d ** (j - i) - X)),
                                             (p * max_pay_off_list[i + 1] + (1 - p) * max_pay_off_list[i]) * df)
print("\n")
print("Option Value from CRR is:",max_pay_off_list[0])    
print("\n","But have to apply interpolation method from Vellekoop")