For J and L (J<L) I want to find that values of p and q which satisfy all the 7 conditions (mention in code) for each 0<=a<b<c<=J and 0<=d<e<f<=L

J=input('J: ');
L=input('L: ');
[a,b,c,d]=ndgrid(1:J,1:J,  1:L, 1:L);
  k=a<b & c<d;
[a,b,c,d]=deal( a(k),  b(k), c(k), d(k)); %all allowed combinations
p_q=cell(6,6);
for j = L:3*L
    for i = 2:j-1
        k =  all( mod( (i.^b-i.^a).*(i.^d-i.^c) ,j)  )        
        p_q{i,j} =[i(k), j(k)]     
    end
end
p_q=vertcat(p_q{:})

The above code is working for single equation but now the condtion is differrent. Here, we have 0<a<b<c<=J, 0<d<e<f<=L;

j=L : uptil the required value is not found

i=2 : j-1

and the conditions are:

for the first equation we have to check for only 0<a<b<=J, 0<d<=L.

mod( (i.^b-i.^a).*(i.^d-i.^c) ,j)

but for the remaining 6 conditions we have to check for 0<a<b<c<=J, 0<d<e<f<=L

and the conditions are

1. mod(((i.^d).*(i.^a-i.^c)+(i.^e).*(i.^c-i.^b)+(i.^f).*(i.^b-i.^a)),j )~=0
2. mod(((i.^d).*(i.^a-i.^c)+(i.^f).*(i.^c-i.^b)+(i.^e).*(i.^b-i.^a)),j )~=0
3. mod(((i.^d).*(i.^b-i.^c)+(i.^f).*(i.^c-i.^a)+(i.^e).*(i.^a-i.^b)),j )~=0
4. mod(((i.^d).*(i.^a-i.^b)+(i.^e).*(i.^b-i.^c)+(i.^f).*(i.^c-i.^a)),j )~=0
5. mod(((i.^d).*(i.^a-i.^b)+(i.^f).*(i.^b-i.^c)+(i.^e).*(i.^c-i.^a)),j )~=0
6. mod(((i.^d).*(i.^b-i.^c)+(i.^e).*(i.^c-i.^a)+(i.^f).*(i.^a-i.^b)),j )~=0

for example

for J=3 and L=4 if we take p=2 q=9 then 6th condition is not satisfied for a=1 b=2 c=3 and d=2 ,e=2, f=3 so we cannot select these values of p and q

next we take p=3 and q=9 then all conditions are not satisfies again we go to next one

on contiuing in this way we will check for all values of p=2 to 9 for q=9 but all the combinaitons are rejected.

and then we will check for q=10, 11, ....... and so on uptil we get the required values of p and q.

Kindly help me in this matter.