...

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!