npts = 7;
deg = 3;
ys = Partition[RandomReal[{0, 1}, npts], deg, deg - 1];
xs = Partition[Sort@RandomReal[{0, 1}, npts], deg, deg - 1];
lxs = Length@xs;
w = Array[x^#2 p@## &, {lxs, deg + 1}, {1, 0}];
t = Thread[Total /@ w == y];
m = Map[{x -> #[[1]], y -> #[[2]]} &,
Transpose /@ Transpose@{xs, ys}, {2}];
joints = (x -> #) & /@ First /@ Rest@xs;
d0 = Equal[#, 0] & /@ D[Differences[Total /@ w], x];
s = Solve[(a = (And @@ Flatten[
ReplaceAll[#[[1]], #[[2]]] & /@ Transpose[{t, m}], 1] &&
And @@ Flatten@Thread[ReplaceAll[d0, joints]] &&