p[x_, y_] := Module[{x1, y1, n, i},
n = Length[x] - 1;
x1 = x;
y1 = y;
h = b = Table[0, {n}];
m = Table[0, {n + 1}];
u = v = Table[0, {n - 1}];
s = Table[0, {n}, 4];
h[[1]] = x1[[2]] - x1[[1]];
b[[1]] = (y1[[2]] - y1[[1]])/h[[1]];
For[i = 2 , i ≤ n, i++,
h [[i]] = x1[[i + 1]] - x1[[i]];
b[[i]] = (y1[[i + 1]] - y1[[i]])/h[[i]], {i, 2, n - 2};]
u[[1]] = 2*(h[[1]] + h[[2]]);
v[[1]] = 6*(b[[2]] - b[[1]]);
For[i = 2 , i ≤ n, i++,
u[[i]] = 2*(h[[i]] + h[[i + 1]]);
v[[i]] = 6*(b[[i + 1]] - b[[i]]);];
m[[1]] = 0;
m[[n + 1]] = 0;
For[ i = n, 2 ≤ i, i--,
m[[i]] = (v[[i]] - h[[i + 1]]*m[[i + 1]])/u[[i]];]]
Coef := Module[{i},
For[ i = 1, i ≤ n, i++,
s[[i, 1]] = y1[[i]];
s[[i, 2]] = b[[i]] - (h [[i]] *(2 m[[k]] + m[[k + 1]]))/6;
s[[i, 3]] = m[[i]]/2;
s[[i, 4]] = (m[[i + 1]] - m[[i]])/6*h[[k]]; ]; ];
... é a vidinha!