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

You can use the inert form  Diff  that will not be calculated, but can be calculated at any time with the command  value :

u[1](T[1], T[2], T[3]):=exp(T[1]*I):
Diff(diff(u[1](T[1], T[2], T[3]), T[1]), T[2]);
value(%);

Of course, you can make inert differentiation for both variables:

u[1](T[1], T[2], T[3]):=exp(T[1]*I):
Diff(u[1](T[1], T[2], T[3]), T[1], T[2]);
value(%);

Here are 2 customised legends with vertical alignment (text rotation is available only from Maple 2018):

restart;
with(plots):
Plots:=plot([t^2,t^3], t=-2..2, color=[blue,red], linestyle=[1,3], view=-2..2, labels=[t,"x,y"]):
Lines:=plot([[1.6,t,t=0.9..1.3],[1.9,t,t=0.9..1.3]], color=[blue,red], linestyle=[1,3]):
Labels:=textplot([[1.55,1.55,x=t^2,rotation=Pi/2],[1.85,1.55,y=t^3,rotation=Pi/2]]):
display(Plots, Lines, Labels);

Use  seq  command for this.

Example:

seq(sin(x), x=0..evalf(Pi/4), 0.1);

     0, 0.9983341665e-1, .1986693308, .2955202067, .3894183423, .4794255386, .5646424734, .6442176872

Check:=proc(n)
local L;
L:=ifactors(n);
if n>1 and type(n,odd) and (not isprime(n)) and `and`(map(t->is(t[2]=1),L[2])[]) and `or`(map(t->is(irem(n-1,t[1]-1)=0),L[2])[]) then true else false fi;
end proc: 

Example of use:
Check(561);
                             true

The previous answer skipped checking of square free.

restart: 
with(plots):
animate(plot, [[sqrt(x)-1, sqrt(x)], x = 0..t, filled = true, color=["Blue", "Red"], view = [0..20, 0..5]], t = 0..20);

Edit.

with(plots):
r:=<cos(t), sin(t), t>:
dr:=diff(r, t):
A:=spacecurve(r, t = -5 .. 5, color=black):
B:=seq(spacecurve(r+s*dr, s = -2 .. 2, color=black), t = -5 .. 5, 0.1):
display(A,B, view=[-1..1,-1..1,-5..5]);

Addition. We can also animate this:

with(plots):
r:=<cos(t), sin(t), t>:
dr:=diff(r, t):
A:=spacecurve(r, t = -5 .. 5, color=black):
B:=b->seq(spacecurve(r+s*dr, s = -2 .. 2, color=black, thickness=0), t = -5 .. b, 0.1):
animate(display,['B'(b)], b=-5..5, frames=100, background=A, view=[-1..1,-1..1,-5..5], orientation=[60,55]);

P:= proc(n::posint)
option remember;
   if n=1 then return 1 fi;
   if type(P(n-1)+1,prime) then return P(n-1)+1 else
   if type(P(n-1),even) then return (n-1)+2 else P(n-1)+3 fi; fi;
end proc:

seq(P(n), n=1..20);
    1, 2, 3, 6, 7, 10, 11, 14, 10, 11, 14, 13, 16, 17, 20, 17, 20, 19, 22, 23

The code above exactly corresponds to the conditions set by you in the first post.

To get the sequence you suggest, you must write

P:= proc(n::posint)
option remember;
   if n=1 then return 1 fi;
   if type(P(n-1)+1,prime) then return P(n-1)+1 else
   if type(P(n-1),even) then return P(n-1)+2 else P(n-1)+3 fi; fi;
end proc:

seq(P(n), n=1..20);
         1, 2, 3, 6, 7, 10, 11, 14, 16, 17, 20, 22, 23, 26, 28, 29, 32, 34, 36, 37

Edit.

Eq:=diff(y(x),x)=cos(x):
ics:=y(0)=0:
Sol:=dsolve({Eq, ics}, numeric):
P:=proc(s)
plot(eval([[x,y(x)]], Sol(s)), style=point, symbol=solidcircle, symbolsize=15);
end proc:
plots:-animate(plots:-display,['P'(s)], s=0..Pi, scaling=constrained, size=[800,300]);

restart;
Expr:=(2*I*Zeta*omega*omega[0] - omega^2 + omega[0]^2)*fourier(omega[r](t), t, omega) = (omega*B[1]*I + B[0])*fourier(delta(t), t, omega);
Expr/op([2,2],Expr)/op([1,1],Expr);

The procedure  Check  tests this hypothesis for the four numbers starting with  n .

Download Conjecture.mw