I suspect that when viewed "over an arbitrary field" (as you say in your Question), the 16x16 case is intractable. But I also suspect that when these matrices are used in practice, the field is small and of characteristic 2. This substantially simplifies the computation.

@WarmtenetwerkBEP All of the other Answers (Kitonum, Joe, VV, Acer) will return true if some vertex occurs in some edge, and false otherwise. Apparently, no-one else had difficulty interpreting the quantifications in your Question.

@Rariusz Does V vary with time, t? If so, then it must be considered as a function by dsolve even if it's only specified by numeric data (that's where interpolation comes in).

If, on the other hand, each record of the file (a record being a group of the four numeric values) is specifying a final value of t where the voltage has been constant from time 0 to time t, then you may consider V as a parameter.

@Rariusz It would be a mistake for you to "close this discussion". Once you specify it precisely, we will be able to solve this problem (or perhaps say that more numeric data is required). You'll not find better help anywhere else.

You need to clarify your requirements using more precise quantifiers, which are words or phrases such as those in my title. Your use of the word "one" is especially confusing because it could mean either "any" or "exactly one". So, do you want to determine whether

1. every vertex occurs in some edge?
2. every vertex occurs in exactly one edge?
3. some vertex occurs in some edge?
4. some vertex occurs in exactly one edge?
5. something else?

@vv Yes, seq is usually faster than an equivalent construction with map and/or elementwise operators. This saddens me, because seq's syntax is ugly. For example, seq requires an index variable (in this context).

@Kitonum There's no need for the other return either.

@Christian Wolinski This degree concept is related to exponents, of course, but it's not about the analytic exp or ln functions; it's a purely algebraic concept. And this is not specific to Maple. If you closely read this Wikipedia article on Euclidean domains (and  "read between the lines") you should see that while degree(0) = 0 would be allowed by the definition, it would be inconvenient, because to effectively perform division-with-remainder we want a function deg such that given any nonzero a and b in the domain there exists q and r in the domain such that  a = q*b + r and deg(r) < deg(b). This ensures termination of Euclid's GCD algorithm, among other things. So, if you think about it for a while, you'll see that we need deg(0) < 0. So, deg(0) = -1 would work, but for reasons already given by VV and me, deg(0) = -infinity is more convenient.

@vv   AFAIK, numeric dsolve never automatically chooses the method other than distinguishing among IVPs, BVPs, and DAE IVPs. Do you have evidence otherwise?

If it's an IVP that hasn't explicity been declared stiff, and no method is declared, then the method will be rkf45.

@ajfriedlan   If this rkf45 solution runs too slow for you (due to the high accuracy requirement), it may be worthwhile specifying another method, such as dverk78.

@Subhan331 The current example is a nonlinear function. Try changing the line where the function is defined. 

And why would you want to directly specify the number of iterations? I could see wanting to limit the maximum number of iterations to prevent infinite loops, but what's the point of extra iterations after it's converged to your desired error tolerance? 

@Subhan331 Is there something wrong with the code that you posted above? I could suggest some minor improvements, but it is a correct implementation of Newton's method.

@mehdibaghaee The ^+ is a transpose operator which I guess doesn't work in your version of Maple or in 2D Input. You can replace it with ^%T

plot([seq(PPP[[1,k],..]^%T, k= 1..rhs(rtable_dims(PPP)[1]))], color= black);

@mehdibaghaee Redundant or tedious commands are almost never needed in Maple. Your sequence of commands above can be reduced to the single line

plot([seq(PPP[[1,k],..]^+, k= 2..rhs(rtable_dims(PPP)[1]))], color= black);

Note that there's no need to specify with(plots), the number of plots (24 in this case), or the number of points per plot (g in this case).

@phbmc What range of values of b are you interested in? For b < 10^5 (approximately), the "brute force" method works well, and I'll post some code soon. That just means trying all possible a1 and checking if the corresponding a2 is an integer in the correct range.

@Subhan331 The code says that the loop ends when the absolute error is less than tol. It just so happens in this case that that happens on the sixth iteration. If you want more iterations, use a smaller value of tol.

As -log[10](tol) gets close to Digits, you will also need to increase Digits to continue to have meaningful iterations. Newton's method tends to double the number of correct Digits on each iteration once it's zeroed in on a target.