It seems the method _d01ajc is only suitable for real-valued functions. Remove this option and the error disappears. In your example, some logarithms will be from negative numbers. Below is a simple example with a similar error:

evalf(Int(ln(x-2), x=0..1, method=_d01ajc));
evalf(Int(ln(x-2), x=0..1));

Error, (in evalf/int) unable to obtain a real result
                  0.3862943611 + 3.141592654 I
 

We have to give names for each root. Replace the last line of your code with the following line:

assign(seq(z[i] = sol1[i], i = 1 .. 4));

You can do this programmatically as in the example below:

restart;
Sys:={x+y-z = 3, x-y-z = 5, -x-y-z = 7};
%piecewise(``, Sys[1], ``, Sys[2], ``, Sys[3]);

In the example below, we first generate a list of 30 random numbers, each in the range from 0 to 5, and then select a sublist of numbers from 0 to 1 from it, indicating the indices:
 

restart;
r:=rand(0...5.):
y:=[seq(r(), n=1..30)];
[seq(`if`(y[n]>0 and y[n]<1,[n,y[n]],NULL), n=1..nops(y))];

  y := [0.2499109628, 2.078923730, 2.252410250, 0.3475034102, 

    1.766288344, 3.629756492, 3.635102067, 3.892007975, 

    3.320931315, 0.008823971507, 3.465092938, 1.209628080, 

    1.722913437, 3.215511858, 4.467963158, 0.5108327385, 

    0.2750667362, 3.363685982, 2.075595849, 1.224136121, 

    3.783826445, 0.5091763968, 3.694547905, 0.2999657008, 

    2.685906488, 1.617331965, 0.2915647810, 1.493696599, 

    4.158416643, 3.518466880]
 [[1, 0.2499109628], [4, 0.3475034102], [10, 0.008823971507], 

   [16, 0.5108327385], [17, 0.2750667362], [22, 0.5091763968], 

   [24, 0.2999657008], [27, 0.2915647810]]
 

restart;
with(Statistics):
x := RandomVariable(Binomial(45, 0.9)):
for i from 0 to 45 do
    P[i]:=[i,ProbabilityFunction(x, i)];
end do:
Points:=convert(P, list);
plot(Points, style=pointline, size=[800,400], labels=["i","x"],labelfont=[times,16]);

   
Edit.

You don't need the combinat package for this. You can just use a nested loop for this (triple loop for your example since k = 3):

restart;
L:=[1,2,3,4,a]: n:=nops(L):
for i1 from 1 to n-2 do
for i2 from i1+1 to n-1 do
for i3 from i2+1 to n do
print([L[i1],L[i2],L[i3]]);
od: od: od:

                           [1, 2, 3]
                           [1, 2, 4]
                           [1, 2, a]
                           [1, 3, 4]
                           [1, 3, a]
                           [1, 4, a]
                           [2, 3, 4]
                           [2, 3, a]
                           [2, 4, a]
                           [3, 4, a]
 

Try the following. In Maple 2018, both methods lead to the same result:

restart;
f := x*y:
JH := int(Heaviside(f-1/2), [x=0..1, y=0..1]); 
JP := int(piecewise(x*y>1/2, 1, 0), [x=0..1, y=0..1]);

restart;
y1:=2*x-4: y2:=-2*x+4:
A:=plot([y1,y2], x=0..4, -2..2, color=blue, thickness=2, labels=[``,``]):
B:=plottools:-sector([2, 0], 0.5, arctan(2) .. Pi-arctan(2), color = pink):
T:=plots:-textplot([[4.3,-0.2,x],[-0.2,2.3,y]], font=[times,16]):
plots:-display(A,B,T, view=[-0.2..4.3,-2..2.3], scaling=constrained);

Or

min(op(xLp_tmp));

restart;
simplify(-0.27059805007310 * sin(0.12 + epsilon) +0.27059805007310 * sqrt(1.-cos(0.12 + epsilon) ^ 2)) assuming sin(0.12 + epsilon)>=0;

                                                           0.

or

expr:=-0.27059805007310 * sin(0.12 + epsilon) +0.27059805007310 * sqrt(1.-cos(0.12 + epsilon)^2);
subs(0.12 + epsilon=t, expr);
simplify(%) assuming t>=0 and t<=Pi;

restart;
G:=GraphTheory:-RandomGraphs:-RandomGraph([seq(v[i],i=1..8)], 10, connected);
GraphTheory:-DrawGraph(G);

Edited.

Using nested loops, we find 1170 solutions. Your example is on the 578th place:

restart;
f:=x->(a*x^2 + b*x + c)/(d*x^2 + e*x + m):
g:=x->d*x^2 + e*x + m:
k:=0:
for a from 1 to 10 do
for b from 1 to 10 do
for c from 1 to 10 do
for d from 1 to 10 do
for e from 1 to 10 do
for m from 1 to 10 do
N:=numer(diff(f(x),x)):
A:=coeff(N,x^2); B:=coeff(N,x); C:=coeff(N,x,0);
disc:=B^2-4*A*C:
if A<>0 and disc>0 and type(sqrt(disc),integer) then x1:=(-B+sqrt(disc))/2/A: x2:=(-B-sqrt(disc))/2/A;
if type(x1,integer) and type(x2,integer) and g(x1)<>0 and g(x2)<>0 then if type(f(x1),integer) and type(f(x2),integer)   then
k:=k+1; L[k]:=[a,b,c,d,e,m] fi; fi; fi;
od: od: od: od: od: od:
L:=convert(L, list):
nops(L);
ListTools:-Search([5,8,2,2,6,5],L);  
                            

                                                                        1170
                                                                         578

restart;
z1:=x1+I*y1: z2:=x2+I*y2: z3:=x3+I*y3:
simplify(evalc(abs(z1+z2+z3)), evalc(map(t->t^2,{abs(z1)=1,abs(z2)=2,abs(z3)=3,abs(9*z1*z2+4*z1*z3+z2*z3)=12}))):
simplify(%);

                                                               2

These 2 files are not loaded. See below my program that solves the problem in a couple of seconds:

restart;
p:=2: k:=0:
while k<10000 do
if isprime(p+36) then k:=k+1; L[k]:=p fi;
p:=nextprime(p);
od:
L:=convert(L,list):

nops(L);
L[-1];

"Warning" is not a bug, but a warning. Your real mistake is that it should be  L[k*n]  and  L[n]  instead of  L__k*n  and  L__n. This is not visible in 2d math and I had to convert your code to 1d math to find and fix this error of yours.

restart; 
Eratosthenes := proc(N::posint) 
local L, primeslist, n, k; 
description "Calculate all primes less than or equal to N"; 
L := Array(2 .. N, i->true); 
for n from 2 to trunc(sqrt(N)) do 
if L[n] = true then 
for k from n while k*n <= N do 
L[k*n] := false od; fi; 
od; 
primeslist := NULL; 
for n from 2 to N do 
if L[n] = true then 
primeslist := primeslist, n fi; 
od; 
primeslist;
end proc:

Eratosthenes(32);

                             2, 3, 5, 7, 11, 13, 17, 19, 23, 29, 31