@Alfred_F 
prove by induction ....

@Math-dashti 

restart;
with(DEtools);
with(plots);
f := (x, y) -> 2*(x - 1)*sin(Pi*y) + (x - 1)^2;
g := (x, y) -> y^2 - y;
lin_sys1 := [diff(u(t), t) = 0, diff(v(t), t) = -v(t)];
lin_sys2 := [diff(u(t), t) = 0, diff(v(t), t) = v(t)];
orig_sys := [diff(x(t), t) = f(x(t), y(t)), diff(y(t), t) = g(x(t), y(t))];
plot_orig1 := DEplot(orig_sys, [x(t), y(t)], t = 0 .. 5, x = 0.5 .. 1.5, y = -0.5 .. 0.5, arrows = SLIM, dirfield = [20, 20], title = "Nonlinear system near (1, 0)", axes = boxed);
plot_orig2 := DEplot(orig_sys, [x(t), y(t)], t = 0 .. 5, x = 0.5 .. 1.5, y = 0.5 .. 1.5, arrows = SLIM, dirfield = [20, 20], title = "Nonlinear system near (1, 1)", axes = boxed);
plot_lin1 := DEplot(lin_sys1, [u(t), v(t)], t = 0 .. 5, u = -0.5 .. 0.5, v = -0.5 .. 0.5, arrows = SLIM, dirfield = [20, 20], title = "Linearized system at (1, 0) → origin (u,v)", axes = boxed);
plot_lin2 := DEplot(lin_sys2, [u(t), v(t)], t = 0 .. 5, u = -0.5 .. 0.5, v = -0.5 .. 0.5, arrows = SLIM, dirfield = [20, 20], title = "Linearized system at (1, 1) → origin (u,v)", axes = boxed);
display([plot_orig1, plot_lin1, plot_orig2, plot_lin2], title = "Comparison: Nonlinear vs. Linearized Systems at Critical Points", insequence = false, layout = [[1, 2], [3, 4]]);


with(LinearAlgebra);
with(DEtools);
with(plots);
fx := (x, y) -> 2*(x - 1)*sin(Pi*y) + (x - 1)^2;
fy := (x, y) -> y^2 - y;
solve({fx(x, y) = 0, fy(x, y) = 0}, {x, y});
J := Matrix(2, 2, [diff(fx(x, y), x), diff(fx(x, y), y), diff(fy(x, y), x), diff(fy(x, y), y)]);
J10 := eval(J, [x = 1, y = 0]);
J11 := eval(J, [x = 1, y = 1]);
print("Jacobiaan bij (1,0):", J10);
print("Jacobiaan bij (1,1):", J11);
sys_uv_10 := {diff(u(t), t) = Pi*v(t), diff(v(t), t) = -v(t)};
sys_uv_11 := {diff(u(t), t) = -Pi*v(t), diff(v(t), t) = v(t)};
field_orig_10 := DEplot([diff(x(t), t) = fx(x(t), y(t)), diff(y(t), t) = fy(x(t), y(t))], [x(t), y(t)], t = 0 .. 5, x = 0.5 .. 1.5, y = -0.5 .. 0.5, dirfield = [20, 20], arrows = SLIM, title = "Niet-lineair systeem rond (1,0)");
field_orig_11 := DEplot([diff(x(t), t) = fx(x(t), y(t)), diff(y(t), t) = fy(x(t), y(t))], [x(t), y(t)], t = 0 .. 5, x = 0.5 .. 1.5, y = 0.5 .. 1.5, dirfield = [20, 20], arrows = SLIM, title = "Niet-lineair systeem rond (1,1)");
field_lin_10 := DEplot([op(sys_uv_10)], [u(t), v(t)], t = 0 .. 5, u = -0.5 .. 0.5, v = -0.5 .. 0.5, dirfield = [20, 20], arrows = SLIM, title = "Linearisatie rond (1,0): u' = πv, v' = -v");
field_lin_11 := DEplot([op(sys_uv_11)], [u(t), v(t)], t = 0 .. 5, u = -0.5 .. 0.5, v = -0.5 .. 0.5, dirfield = [20, 20], arrows = SLIM, title = "Linearisatie rond (1,1): u' = -πv, v' = v");
display([field_orig_10, field_lin_10, field_orig_11, field_lin_11], title = "Origineel versus lineair systeem rond kritieke punten");

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