Mat[r_,rank_]:=(For[a={};b={};c=Length[r];d=0,d<c,b=Union[b,Subsets[r[[d]],{rank}]];
a=Union[a,Subsets[r[[d]],{rank-1}]],d=d+1];
Table[Table[If[Intersection[b[[i]],a[[j]]]==a[[j]],-(-1)^(Position[b[[i]],Complement[b[[i]],a[[j]]][[1]]][[1]])[[1]],0],{i,1,Length[b]}],{j,1,Length[a]}])

soln[r_,rank_,q_]:=(z=Mat[r,rank];
Table[m[i],{i,1,Length[b]}]/.Solve[z.Transpose[{Table[m[i],{i,1,Length[b]}]}]==Transpose[{Table[0,{i,1,Length[a]}]}],Table[m[i],{i,1,Length[b]}],Modulus->q][[1]])

base[r_,rank_,q_]:=(w=soln[r,rank,q];j=1;
While[(w/.Table[C[i]->Switch[i==j,False,0,True,1],{i,1,Length[b]+1}])≠Table[0,{k,1,Length[b]}],j=j+1];
j-1)

basis[r_,rank_,q_]:=(p=soln[r,rank,q];ba=base[r,rank,q];
Table[p/.(Table[C[i]->Switch[i==j,False,0,True,1],{i,1,ba}]),{j,1,ba}])

doesithave[r_,rank_,q_]:=(If[base[r,rank,q]==0,False,h=0;
ba=base[r,rank,q];
bas=basis[r,rank,q];A={};m=True;
While[m&&Length[A]≤q^ba/(q-1),For[s=1;e={},s≤ba,s=s+1,e=Append[e,{Random[Integer,{0,q-1}]}]];
If[Not[MemberQ[A,Transpose[e][[1]]]],f=Mod[Transpose[bas].e,q];A=Append[A,Transpose[e][[1]]]];
If[Product[f[[i]],{i,1,Length[f]}][[1]]≠0,h=1;Break[]]];
If[h==1,True,False]])

doesithave[{{1,2},{2,3},{3,4},{4,5},{5,1},{1,6},{2,7},{3,8},{4,9},{5,10},{6,8},{7,9},{8,10},{9,6},{10,7}},2,4] (*This is the 4-flows on the Petersen Graph problem*)

doesithave[{{1,2},{2,3},{3,4},{4,5},{5,1},{1,6},{2,7},{3,8},{4,9},{5,10},{6,8},{7,9},{8,10},{9,6},{10,7}},2,5](*This is the 5-flows on the Petersen Graph problem*)

doesithave[{{1,2,3,4,5,6}},4,6](*This is the simplex Brad was talking about - it has 6 flows (above), but it has no 5 flows (below).*)

doesithave[{{1,2,3,4,5,6}},4,5]

(*The Mat function takes the simplicial complex, entered as lists of "faces", and the rank of the higher face (2 in the case of a normal graph, 3 in the case of the tetrahedron, and so on), and outputs the matrix to go with it.*)

(*The soln function takes in the complex, the rank and the value for q, and basically row-reduces the matrix, computing the form of a general solution.*)

(*The base function just calculates the number of free variables of the general solution.*)

(*The basis function takes each free variable in the general solution and sets it equal to 1 and all others equal to 0, and tabulates each of these solutions as a matrix.  Basically this is just a roundabout way of calculating the basis of solutions so that we can take linear combinations of them mod q and see whether they're nonzero.*)

(*The doesithave function again takes in the complex, rank and q and determines if it has a q flow in the rank'th dimension.  It does this nondeterministically, by randomly choosing a tuple of the correct size and with integers between 0 and q-1 and seeing if multiplying this by the basis matrix gives a nonzero vector.  If so, it breaks automatically and returns true.  If not, it adds the randomly generated vector to a stockpile called A and continues on, at each step checking if A contains the randomly generated vector.  It stops with a "NO" when the size of A is at least 1/(q-1) of the total number of possible vectors. The speed of searching through A is not really a hindrance, so if we wanted a sure-fire program I could change it so that it starts with the set of all possible vectors and randomly picks one, tests it, and removes it.*)