We can easily solve this problem using the  seq  command with the  `if` - condition:

L:=[$ 1..50];
L1:=[seq(`if`(irem(i,2)=0,NULL, L[i]), i=1..nops(L))];
L2:=[seq(`if`(irem(i,3)=0,NULL, L1[i]), i=1..nops(L1))];

Edit. This method is about 2 times faster than using the  remove  command:

L:=[$ 1..1000000]:
CodeTools:-Usage([seq(`if`(irem(i,2)=0,NULL, L[i]), i=1..nops(L))]):
CodeTools:-Usage(RemoveEveryNth(L, 2)):

I have simplified your procedures somewhat by removing unused formal parameters and unnecessary local variables. Note also that the external call of  packages does not work in the procedures. Use  the  uses  command to call the package inside the procedure.

Download func_new.mw

1. First, you made a common syntax error: the exponent is encoded in Maple as  exp(x)  not  e^x  (in Maple  e  is only a symbol rather than the number  2.71828...).

2. Your integral contains 4 parameters (a, b, c, d), but their role is different. The most important is the parameter  b and we can easily get a solution (in many cases) depending on the value of  this parameter. The examples below show that the solution form depends heavily on the value of  b :

restart;
f:=b->d*exp(-b*x)/(((a*exp(-2*b*x)+c*exp(-4*x))));
F:=b->int(f(b),x);

 # Examples:
for b from 0 to 5 do
F(b);
od; 

Edit.

plots:-implicitplot((x-2)^2*(x-3)+epsilon^2=0, x=-1..6, epsilon=-5..5, gridrefine=4);
solve((x-2)^2*(x-3)+epsilon^2=0, x, parametric, real);

So we have 3 values for epsilon:  -2*sqrt(3)*(1/9),  0,  2*sqrt(3)*(1/9)

Use the  eliminate  command for this. The reverse representation is ambiguous and depends on the choice of parameters. For your example, if we take  x  and  y  as parameters, it is trivial:

Download elim.mw

A:=<1,2; 3,4>;
B:=<5,6; 7,8; 9,10; 11,12>;
<A[1], B[2], B[4]>;

An alternative to the  map  command would be the  seq  command:

restart;
L:=[[1,2,3],[7,8,9],[13,12,11]]:
[seq(l[3], l=L)];
L1:=[1, [2, 3], [4, [5, 6], 7], [8, 3], 9]:
[seq(`if`(l::list,l[1],l), l=L1)];

Addition. You do not need to form a system of equations in any form. Using the  LinearAlgebra:-LinearSolve  command you can solve the system simply by specifying the matrix of coefficients in front of the unknowns and the vector of the right side:

M:=(Matrix([[rho*h*omega[m,n]^2+k[1, 1],k[1,2],k[1,3],k[1,4]],[k[2,1],rho*h*omega[m,n]^2+k[2,2],k[2,3],k[2,4]],[k[3,1],k[3,2],k[3,3]+rho*h*omega[m,n]^2,k[3,4]],[k[4,1],k[4,2],k[4,3],k[4,4]+rho*h*omega[m,n]^2]]));
V:=-E*(Vector(1..4,[k[1,5],k[2, 5],k[3,5],k[4,5]]));
LinearAlgebra:-LinearSolve(M, V);

Download solve_(1)_new.mw

The  AlmostCarmichaelNumbers  procedure finds all the numbers that meet your requirements for a number  k  and for a range  R  (which you specify when using) .

AlmostCarmichaelNumbers:=proc(k::posint,R::range)
local T, N, m, n, L;
T:=table(); N:=0; m:=lhs(R);
for n from `if`(m::odd,m,m+1) by 2 to rhs(R) do
L:=ifactors(n);
if nops(L[2])=k then if `and`(map(t->is(t[2]=1),L[2])[]) and `+`(seq(`if`(irem(n-1,t[1]-1)=0,1,0),t=L[2]))=k-1 then N:=N+1; T[N]:=n fi; fi;
od;
convert(T, list);
end proc:

Examples of use for k=3, 4, 5 :

AlmostCarmichaelNumbers(3,1..10000);
AlmostCarmichaelNumbers(4,1..10000);
AlmostCarmichaelNumbers(5,1..10000);
              

A:=(-108*q - 12*sqrt(12*p^3 + 81*q^2))^(1/3)/6;
B:=(A^3)^(1/3);
c:=op([1,2], B);
subsop([1,2]=-(c^2)^(1/2), B);

Addition. This method can be used in some cases when simplifying radicals, for example:

restart;
A:=sqrt(11-4*sqrt(6))+sqrt(21+6*sqrt(6));
sqrt(A^2);

The differential equation must be entered with explicit indication of the derivative (not differentials) using Maple syntax. For your example, you can do this:

restart;
de:=2*y*dx+3*x*dy=0;
subs([dy=diff(y(x),x)*dx,y=y(x)], de);
dsolve(%, y(x));

Just differentiate your equation 2 times, taking  y as a function of  x :

restart;
Curves:=x^2+A*y(x)^2=B:
Eq:=diff(Curves,x,x);
dsolve({Eq,y(0)=sqrt(B/A),D(y)(0)=0}, y(x), implicit);
solve(%, B)=B;

Edit. Sure enough the right solution is Rouben's one, not me. I forgot to eliminate the constant  A . Rouben very politely "did not notice" this error.

You can do this (like as in the original post) without  DataFrame  , using only matrices:

A:=Matrix(3, (i,j)->i);
DF:=< <<D1|D2|D3>, A> | <``,D1,D2,D3> >;

You can use the  Explore  command for this:

Explore(plot([sqrt(a*x^3+b*x^2+c*x+d),-sqrt(a*x^3+b*x^2+c*x+d)], x=-3..3, color=red, view=-7..7), a=-3...3., b=-3...3., c=-3...3., d=-3...3.);

Maple's possibilities are not limitless, so just reduce the ranges for calculations. It is easy to check that one of the roots is exact 0:

Download calculs_new.mw