Working on cubic spline interpolation, can`t figure out why all my double[] are full of NaN values...

let thrd (_, _, c) = c

let frst (c, _, _) = c

let scnd (_, c, _) = c

let LUdecomp3 (c : double[], d : double[], e : double[]) =

let n : int = d.Length

for k = 1 to n-1 do

let lam = c.[k-1] / d.[k-1]

d.[k] <- d.[k] - lam * e.[k-1]

c.[k-1] <- lam

(c, d, e)

let LUsolve3 (c : double[], d : double[], e : double[], b : double[]) =

let n = d.Length

for k = 1 to n-1 do

b.[k] <- b.[k] - c.[k-1] * e.[k-1]

b.[n-1] <- b.[n-1] / d.[n-1]

for k = n-2 to -1 do

b.[k] <- (b.[k] - e.[k] * b.[k+1]) / d.[k]

b

let curvatures (xData : double[], yData : double[]) : double[] =

let n : int = xData.Length - 2

let mutable c : double[] = Array.zeroCreate(n+1)

let mutable d : double[] = Array.zeroCreate(n+2)

let mutable e : double[] = Array.zeroCreate(n+1)

let mutable k : double[] = Array.zeroCreate(n+2)

for i = 0 to n-1 do

c.[i] <- xData.[i] - xData.[i+1]

for i = 1 to n do

d.[i] <- 2.0 * (xData.[i-1] - xData.[i+1])

e.[i] <- xData.[i] - xData.[i+1]

k.[i] <- 6.0 * (yData.[i-1] - yData.[i+1])

/ (xData.[i-1] - xData.[i])

- 6.0 * (yData.[i] - yData.[i+1])

/ (xData.[i] - xData.[i+1])

let res = LUdecomp3(c, d, e)

k <- LUsolve3(frst(res), scnd(res), thrd(res), k)

k