hamilton_syteem_plotsyteem_mprimes_25-6-2025.mw

Drag to resize the screen and maybe you will see this in de right upper corner?

restart;
with(plots):
f := x -> cos(Pi*sqrt(x^2 + x + 1));
points := [seq([n, f(n)], n = 1 .. 100)]:
continuousPlot := plot(f(x), x = 1 .. 100, color = red, thickness = 2):
discretePlot := pointplot(points, symbol = solidcircle, color = blue):
linePlot := plottools[curve](points, color = green, linestyle = dot, thickness = 3):
display([continuousPlot, discretePlot, linePlot], title = "Gedrag van cos(Pi*sqrt(x^2 + x + 1)) versus f(n)", labels = ["x", "f(x)"], axes = boxed);

I think you can calculate with "polar" vectors in the Vector Calculus package , by setting up the appropriate coordinate system for this through the package.

But graphically you need to calculate the vectors back to Cartesian coordinates to plot them.

Create a small procedure for this then

There is a polar plot command, but that is for functions in polar coordinates.

from example 1 , as a start 

is almost identical with the LAD_complex procedure , see also plot comparison 

i changed the procedure, must be corrected manual

 procedure_complexe_pde_mprimesVERSIONA_11-6-2025.mw

Hopefully you are satisfied,because it is slightly different, but it can be made the same

>  Analysis completed!

Download Can_plot_the_matrix_instead_of_linear_system_of_odemprimes10-6-2025.mw

?
 

DOLOOP_SUBST_9-6-2025MPRIMES.mw

restart;
with(DEtools);
with(plots);
with(LinearAlgebra);
f := y -> y^2*(1 - y^2);
sys := [diff(x(t), t) = -x(t), diff(y(t), t) = f(y(t))];
trajectories := [[x(0) = 1, y(0) = 1.5], [x(0) = -1, y(0) = 0.5], [x(0) = 0.5, y(0) = -1.5]];
phase_plot := DEplot(sys, [x(t), y(t)], t = 0 .. 10, trajectories, x = -2 .. 2, y = -2 .. 2, arrows = medium, dirfield = [20, 20], linecolor = blue, title = "Phase Plane with Equilibrium Points and Trajectories");
equilibrium_points := pointplot([[0, 0], [0, 1], [0, -1]], symbol = solidcircle, color = black, symbolsize = 15);
display([phase_plot, equilibrium_points], axes = normal, scaling = constrained, title = "Phase Plane with Equilibrium Points and Trajectories");
print("\n### STABILITY ANALYSIS ###");
F := [-x, y^2*(-y^2 + 1)];
J := Matrix([[diff(F[1], x), diff(F[1], y)], [diff(F[2], x), diff(F[2], y)]]);
print("\n1. Analysis for (0,0):");
J0 := eval(J, {x = 0, y = 0});
print("Jacobian at (0,0):");
print(J0);
eig0 := Eigenvalues(J0);
print("Eigenvalues:", eig0);
print("Conclusion:");
print("   - Eigenvalues: λ1 = -1, λ2 = 0");
print("   - One zero eigenvalue → linearization incomplete");
print("   - From original equation dy/dt = y^2(1-y^2):");
print("     * For y ≈ 0: dy/dt ≈ y^2 > 0 → y moves AWAY from 0");
print("   - Conclusion: (0,0) is UNSTABLE (non-attracting)");
print("\n2. Analysis for (0,1):");
J1 := eval(J, {x = 0, y = 1});
print("Jacobian at (0,1):");
print(J1);
eig1 := Eigenvalues(J1);
print("Eigenvalues:", eig1);
print("Conclusion:");
print("   - Eigenvalues: λ1 = -1, λ2 = -2");
print("   - Both eigenvalues negative → stable node");
print("   - All nearby trajectories converge to this point");
print("   - Conclusion: (0,1) is STABLE and ATTRACTING");
print("\n3. Analysis for (0,-1):");
Jm1 := eval(J, {x = 0, y = -1});
print("Jacobian at (0,-1):");
print(Jm1);
eigm1 := Eigenvalues(Jm1);
print("Eigenvalues:", eigm1);
print("Conclusion:");
print("   - Eigenvalues: λ1 = -1, λ2 = -2");
print("   - Both eigenvalues negative → stable node");
print("   - All nearby trajectories converge to this point");
print("   - Conclusion: (0,-1) is STABLE and ATTRACTING");
print("\n### SUMMARY OF EQUILIBRIUM POINTS ###");
print("1. (0,0): Unstable (non-hyperbolic point)");
print("2. (0,1): Stable and attracting node");
print("3. (0,-1): Stable and attracting node");
print("\nNote: The system shows bistability with two stable equilibria at (0,1) and (0,-1)");

complexe_ber_mprimes-8-6-2025.mw

Check it , if you can use it , good luck

more friendly version
complexe_ber_mprimesDeel2met_nummering-8-6-2025.mw

@salim-barzani 
 

restart;
with(inttrans);
with(PDEtools);
with(DEtools);
with(Physics);
declare(u(x, t), quiet);
declare(v(x, t), quiet);
undeclare(prime);
pde := u(x, t) + Physics:-`*`(diff(u(x, t), x $ 2), I) + Physics:-`*`(2, I, diff(Physics:-`*`(u(x, t), conjugate(u(x, t))), x), u(x, t)) + Physics:-`*`(Physics:-`^`(u(x, t), 2), Physics:-`^`(conjugate(u(x, t)), 2), I, u(x, t));
pde_linear, pde_nonlinear := selectremove(term -> not has(eval(term, u(x, t) = T*u(x, t))*1/T, T), expand(pde));
B[0] := -Physics:-`*`(I, Physics:-`^`(u[0], 3), Physics:-`^`(conjugate(u[0]), 2));
B1[0] := -Physics:-`*`(2, I, Physics:-`^`(u[0], 2), diff(u[0](x), x));
T[0] := -Physics:-`*`(2, I, u[0], diff(u[0](x), x), conjugate(u[0]));
P[0] := B[0];
Q[0] := B1[0];
R[0] := T[0];
A[0] := P[0] + Q[0] + R[0];
u[0] := Physics:-`*`(beta, exp(Physics:-`*`(I, x)));
u_conj[0] := Physics:-`*`(beta, exp(-Physics:-`*`(I, x)));
A_eval := subs(conjugate(u[0]) = u_conj[0], A[0]);
uxx := diff(u[0](x), x $ 2);
expr := -Physics:-`*`(I, uxx) + A_eval;
u[1] := invlaplace(Physics:-`*`(laplace(expr, t, s), Physics:-`^`(s, -1)), s, t);
print("u[0] =", simplify(u[0]));
print("u[1] =", simplify(u[1]));

@salim-barzani 
splitsing_pde_conjugate_vorm_4-6-2025mprimes.mw

LaplaceAdomian_Decomposition_Method_(LADM)mprimesDEF-4-6-2025.mw

The first odetest as example 

restart;
with(PDEtools);
with(LinearAlgebra);
with(SolveTools);
undeclare(prime);
declare(G(xi));
declare(U(xi));
ode := diff(G(xi), xi) = ln(d)*(A + B*G(xi) + C*G(xi)^2);
G1 := G(xi) = -B/(2*C) + sqrt(4*A*C - B^2)*tan(1/2*sqrt(4*A*C - B^2)*xi)/(2*C);
H := unapply(rhs(G1), xi);
Hscaled := unapply(H(xi*ln(d)), xi);
Gscaled := G(xi) = Hscaled(xi);
(odetest(Gscaled, ode) assuming (-4*A*C + B^2 < 0, 0 < d, d <> 1, C <> 0));

===========================================================

@Alfred_F 

Download how_to_use_ggd-mprimes19-5-2025.mw