OK, so here is a really scetchy solution (up to numerical accuracy). It can compute the requested vectors not only in the case of 5 vectors, but for arbitrary number (set variable point
to the desired value).
Remove["Global`*"]
point = 7;
SeedRandom;
rnd[x_] := RandomReal[{-10, 10}];
mat = {};
Do[mat = Append[
mat, {Subscript[p0, i], sg Subscript[p1, i], sg Subscript[p2, i],
sg Subscript[p3, i]}];, {i, 1, point}];
vec = {};
Do[vec = Append[vec, 1];, {i, 1, point}];
vec1 = Diagonal[(mat /. sg -> 1).(Transpose[mat] /. sg -> -1)];
vec2 = Transpose[mat /. sg -> 1].vec;
vecf = vec1;
Do[vecf = Append[vecf, vec2[[i]]];, {i, 1, Length[vec2]}];
res = vecf - vecf;
now = 4 point - 4 - point;
now1 = now - point + 1;
If[now1 > -1, cnt3 = point - 1;, If[now > -1, cnt3 = now, cnt3 = 0]];
now = now1;
now1 = now - point + 1;
If[now1 > -1, cnt2 = point - 1;, If[now > -1, cnt2 = now, cnt2 = 0]];
now = now1;
now1 = now - point + 1;
If[now1 > -1, cnt1 = point - 1;, If[now > -1, cnt1 = now, cnt1 = 0]];
sol = mat[[All, 1]];
Do[Subscript[p3, i] = rnd[1];, {i, 1, cnt3}];
Do[sol = Append[sol, Subscript[p3, i]], {i, cnt3 + 1, point}];
Do[Subscript[p2, i] = rnd[1];, {i, 1, cnt2}];
Do[sol = Append[sol, Subscript[p2, i]], {i, cnt2 + 1, point}];
Do[Subscript[p1, i] = rnd[1];, {i, 1, cnt1}];
Do[sol = Append[sol, Subscript[p1, i]], {i, cnt1 + 1, point}];
sub = FindInstance[vecf == res, sol, Reals][[1]];
Do[Subscript[p0, i] = Subscript[p0, i] /. sub;
Subscript[p1, i] = Subscript[p1, i] /. sub;
Subscript[p2, i] = Subscript[p2, i] /. sub;
Subscript[p3, i] = Subscript[p3, i] /. sub;, {i, 1, point}];
The variable vecf
containes all the constraints. After the solution has been found all the entries of vecf
are properly numerically zero.