Judging by Eq. 12 in the paper cited by Burkardt, it does look like a term is missing:
For k_4 (j = 4) the sum in the bracketed part of functions F and G, should go from i=1 to j - 1 = 3. This does seem to imply the line should read
x4 = x1 + a41 * k1 + a42 * k2 + a43 * k3