@salim-barzani 

restart;
with(plots):

# Definieer de functie voor de 3D plot met parameters en componentkeuze
energie_3d_keuze := proc(comp, m, p, q, r, theta, w, v, t, alpha1, alpha2, alpha3, alpha4)
    local Delta, xi, Knum, Pnum, Hnum, plot_expressie, plot_kleur, plot_titel;
    # Bereken Delta
    Delta := sqrt(4*alpha1*w^2 + (4*theta*alpha3 - 4)*w + 3*alpha4*m^2 + 4*theta^2*alpha2 - 4*r) / 
             (2*sqrt(p^2*alpha1 + p*q*alpha3 + q^2*alpha2));
    # Definieer xi
    xi := p*x + q*y - t*v;
    # Bereken energie componenten
    Knum := - (m^2/(16*p^2*alpha1 + 16*p*q*alpha3 + 16*q^2*alpha2)) * 
            (3*m^2*alpha4 + 4*theta^2*alpha2 + 4*theta*w*alpha3 + 4*w^2*alpha1 - 4*r - 4*w) * 
            (-1 + cos(2*Delta*xi));
     Pnum := (1/(4*p^2*alpha1 + 4*p*q*alpha3 + 4*q^2*alpha2)) * cos(Delta*xi)^2 * 
            (alpha4*m^2*cos(Delta*xi)^2 + 2*alpha1*w^2 + (2*theta*alpha3 - 2)*w + 2*theta^2*alpha2 - 2*r) * m^2;
     Hnum := (m^2/(8*p^2*alpha1 + 8*p*q*alpha3 + 8*q^2*alpha2)) * 
            (2*alpha4*m^2*cos(Delta*xi)^4 - 3*alpha4*m^2*cos(Delta*xi)^2 + 
             3*alpha4*m^2 + 4*theta*w*alpha3 + 4*alpha1*w^2 + 4*theta^2*alpha2 - 4*r - 4*w);
    
    # Kies welke component geplot wordt
    if comp = 1 then
        plot_expressie := Knum;
        plot_kleur := red;
        plot_titel := "Kinetische Energie (Knum)";
    elif comp = 2 then
        plot_expressie := Pnum;
        plot_kleur := blue;
        plot_titel := "Potentiële Energie (Pnum)";
    else
        plot_expressie := Hnum;
        plot_kleur := green;
        plot_titel := "Totale Energie (Hnum)";
    end if;
    
    # Maak de 3D plot
    plot3d(plot_expressie, x = -5..5, y = -5..5, 
           view = [-5..5, -5..5, 0..0.35], color = plot_kleur,
           title = sprintf("%s (m=%.2f, p=%.2f, q=%.2f, t=%.2f)", plot_titel, m, p, q, t),
           labels = ["x", "y", "energie"],
           style = patchcontour, axes = boxed);
end proc:

# Maak de Explore versie met componentkeuze (zonder controller optie)
Explore(energie_3d_keuze(comp, m, p, q, r, theta, w, v, t, alpha1, alpha2, alpha3, alpha4),      
        parameters = [
            comp = [1, 2, 3],
            m = 0.1..5.0, p = 0.1..5.0, q = 0.1..5.0,
            r = 0.1..5.0, theta = 0.1..5.0, w = 0.1..5.0,
            v = 0.1..5.0, t = 0.1..5.0,
            alpha1 = 0.1..5.0, alpha2 = 0.1..5.0,
            alpha3 = 0.1..5.0, alpha4 = 0.1..5.0
        ],      
        initialvalues = [
            comp = 1, m = 1, p = 1, q = 1, r = 1, theta = 1, w = 1, v = 1, t = 1,
            alpha1 = 1, alpha2 = 1, alpha3 = 1, alpha4 = 1
        ],    
        placement = right);

Maybe this gui ?

Beta cannot be 0 , as a quick test, not symbolic , but numerically 

# Choose random numerical values for the parameters
#params := {beta = 0.5, eta1 = 1.2, eta2 = 0.8, gamma1 = 0.3, gamma2 = -0.4, N = 2};

# Evaluate both expressions with these values
eval(G0, params);
eval(G01, params);

# Check if they are (approximately) equal
evalb(eval(G0 - G01, params) = 0);
# Or use a tolerance
is(abs(eval(G0 - G01, params)) < 1e-10);

Maybe useful?

nanofluid_figures_mprimesDEFA_19-2-2026.mw

t := 5*Pi/9:
L := [-tan(4*Pi/9) + 4*sin(4*Pi/9)];
evalf~(L,20)[];

identify(-1.7320508075688772935);

A more practical approach ,well done..Maple, for recognizing a root form from the decimal form.

oneliner..( a long line ) 

toPrefix := e -> `if`(type(e,atomic), e,[op(0,e), seq(toPrefix(op(i,e)), i=1..nops(e))]):

ex := (a+b*c)^2:
toPrefix(ex);


exUltra :=
subs(a=b+c,
    expand(
        int(
            diff((x+y)^5*sin(exp(x^2+y^2)),x)
            /(1+sum(k*z^k,k=1..4)),
            x=0..1
        )
    )
):

toPrefix(exUltra );

exInsane :=
((f@g)@@3)(a+b*c)
+
(exp@@2)(sin@@3(x+y))
+
((a&*b)&+c)^5;

 exInsane := f(g(f(g(f(g(b c + a)))))) + @@(exp, 2)(@@(sin, 3))

                   5
    + (a &* b) &+ c 

toPrefix(exInsane );

A name is  a variable when there is a valued assigned to it. 

Used seq command here instead  a loop construction
Note: how fast the trees are growing :-) 

>      >  restart; with(Fractals[LSystem]): with(plots): # Define the L-system generator LGenerator := proc(n, Init, Angle0, Param, Rules, Color)     local Iteration;          # Use try-catch to avoid errors     try         if n = 0 then             Iteration := Init;         else             Iteration := Iterate(Init, Rules, n);         end if;     catch:         # If there's an error, use the initial string         Iteration := Init;     end try;          # Draw the L-system     LSystemPlot(         Iteration,         Param,         initialangle = Angle0,         color = Color,         scaling = constrained,         axes = none,         thickness = 2     ); end proc: # First L-system (number 17) Init17 := "X": Angle17 := 90: Param17 := [   "F" = "draw:1",   "+" = "turn:-25.7",   "-" = "turn:25.7",   "[" = "push:position;push:angle",   "]" = "pop:position;pop:angle" ]: Rules17 := [   "F" = "FF",   "X" = "F[+X][-X]FX" ]: # Second L-system (number 18) Init18 := "X": Angle18 := 90: Param18 := [   "F" = "draw:1",   "+" = "turn:20",   "-" = "turn:-20",   "[" = "push:position;push:angle",   "]" = "pop:position;pop:angle" ]: Rules18 := [   "F" = "FF",   "X" = "F[+X]F[-X]+X" ]: # Test both L-systems without animation print("L-system 17, iteration 4:"); display(LGenerator(4, Init17, Angle17, Param17, Rules17, "Red")); print("L-system 18, iteration 4:"); display(LGenerator(4, Init18, Angle18, Param18, Rules18, "DarkGreen")); # Create smoother animations with more frames print("Smooth animation of L-system 17 (iterations 0-7):"); display([seq(LGenerator(i, Init17, Angle17, Param17, Rules17, "Red"), i = 0 .. 7)],         insequence = true,         frames = 64);  # More frames for smoother animation print("Smooth animation of L-system 18 (iterations 0-7):"); display([seq(LGenerator(i, Init18, Angle18, Param18, Rules18, "DarkGreen"), i = 0 .. 7)],         insequence = true,         frames = 64);  # More frames for smoother animation # Optional: Show all iterations side by side in a grid print("All iterations of L-system 17 (0-7):"); display(Matrix(2, 4, [seq(LGenerator(i, Init17, Angle17, Param17, Rules17, "Red"), i = 0 .. 7)])); print("All iterations of L-system 18 (0-7):"); display(Matrix(2, 4, [seq(LGenerator(i, Init18, Angle18, Param18, Rules18, "DarkGreen"), i = 0 .. 7)])); # Even smoother version with interpolation (more iterations shown) print("Extra smooth animation of L-system 17 (more frames):"); display([seq(LGenerator(i, Init17, Angle17, Param17, Rules17, "Red"), i = 0 .. 7)],         insequence = true,         frames = 128);  # Even more frames for extra smoothness print("Extra smooth animation of L-system 18 (more frames):"); display([seq(LGenerator(i, Init18, Angle18, Param18, Rules18, "DarkGreen"), i = 0 .. 7)],         insequence = true,         frames = 128); >

Download fractal_animatie_mprimes_26-1-2026.mw

@jalal 
Some experiment, but needs more thought 
I do need more math background/information for this figure in order to get this as a plot in Maple 

restart;
with(plots):

# parameters
p := 54.62;           # exponent
N := 120;             # number  rotaties
tmin := -1;
tmax := 1;

curves := []:

for k from 0 to N do
    theta := 2*Pi*k/N;

    x := t*cos(theta) - t^p*sin(theta);
    y := t*sin(theta) + t^p*cos(theta);

    curves := [op(curves),
        plot([x, y, t = tmin .. tmax],
             color = ColorTools:-Color([sin(theta)^2,
                                         cos(theta)^2,
                                         abs(sin(2*theta))]))];
end do:

display(curves,
        scaling = constrained,
        axes = none,
        background = white);

Download chirped_parameter_eta1_eneta_2_mprimes_11-1-2026.mw

@salim-barzani 
gereduceerde_ode_oplosser_mprimes_DEF_11-1-2026.mw

@salim-barzani 

HirotaTau :=
proc(N::posint, xi, A)
    local mu, f0, mus, muvals, f, i, j, mui, muval;

    # binaire variabelen mu[1]..mu[N]
    mu := array(1..N):

    # generieke subset-term
    f0 :=
    exp(
        add(mu[i]*xi[i], i=1..N)
      + add(
            add(mu[i]*mu[j]*ln(A[i,j]),
                i=1..j-1),
            j=1..N
        )
    );

    # alle binaire combinaties van mu (Gray code)
    mus := Iterator:-BinaryGrayCode(N);

    # vervangingen mu[i]=0/1
    muvals :=
    seq(
        { seq(mu[i], i=1..N) =~ seq(mui) },
        mui in mus
    );

    # subset-som
    f := add( eval(f0, muval), muval in muvals );

    return simplify(f);
end proc:

#restart;

N := 3:

xi := table():
A  := table():

f3 := HirotaTau(N, xi, A);

First doing by hand , to see where it goes wrong in the code : 
Step 3 is a vital one , but not easy to get , by knowing this 
A nightmare this 

# Step 1: Define original expression
expr := sqrt(4 - 2*y*(x + 1));

# Step 2: Expand the expression inside the root
expandedjd := expand(4 - 2*y*(x + 1));  # Gives -2 x y - 2 y + 4

# Step 3: Factor using symbolic simplification
factored := simplify(expandedjd, symbolic);  # Returns 2*(2 - y*(x + 1))

# Step 4: Apply the square root to the factored expression
root_factored := sqrt(factored);  # sqrt(2*(2 - y*(x + 1)))

# Step 5: Split the root into a product of square roots
split_root := sqrt(2)*sqrt(2 - y*(x + 1));

# Step 6: Verify equivalence
check := simplify(root_factored - split_root);  # Should return 0

Yes — for elliptic integrals of the second (and first) kind with real arguments, this is always true.
sin(x) never causes branching there; all branch information runs via cos(x). ( ai  ?) 
Turns out it has a very unique approach?
 

I don't know how to use all elliptic integrals to obtain the desired output by specifying conditions for the integral definitions.
That would need to be investigated, but is it now known for one elliptic integral, but not direct?

restart:
f := sqrt(1 + sin(x)^2);
F := int(f, x);
df := diff(F, x):
simplify(df - f);
simplify(df - f, symbolic);

restart;
f := sqrt(1 + sin(x)^2);
F := int(f, x);
F := simplify(F);
#F := simplify(F,symbolic);
df := diff(F, x);
assume(cos(x) > 0):
simplify(df - f);

piecewise(
  cos(x) > 0,  EllipticE(sin(x), I),
  cos(x) < 0, -EllipticE(sin(x), I)
);

restart:
f := sqrt(1 + sin(x)^2):
F := int(f, x);
df := diff(F, x);
simplify(df - f);  # ==> must be give  0