diff --git a/build/pkgs/maxima/package-version.txt b/build/pkgs/maxima/package-version.txt
index b1e9710a45b..41b4e210ed2 100644
--- a/build/pkgs/maxima/package-version.txt
+++ b/build/pkgs/maxima/package-version.txt
@@ -1 +1 @@
-5.38.1
+5.38.1.p1
diff --git a/src/doc/de/tutorial/tour_algebra.rst b/src/doc/de/tutorial/tour_algebra.rst
index bd788d9d7e7..5864cbff79e 100644
--- a/src/doc/de/tutorial/tour_algebra.rst
+++ b/src/doc/de/tutorial/tour_algebra.rst
@@ -212,7 +212,7 @@ Lösung: Berechnen Sie die Laplace-Transformierte der ersten Gleichung
 
     sage: de1 = maxima("2*diff(x(t),t, 2) + 6*x(t) - 2*y(t)")
     sage: lde1 = de1.laplace("t","s"); lde1
-    2*(-%at('diff(x(t),t,1),t=0)+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)
+    2*((-%at('diff(x(t),t,1),t=0))+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)
 
 Das ist schwierig zu lesen, es besagt jedoch, dass
 
@@ -228,7 +228,7 @@ Laplace-Transformierte der zweiten Gleichung:
 
     sage: de2 = maxima("diff(y(t),t, 2) + 2*y(t) - 2*x(t)")
     sage: lde2 = de2.laplace("t","s"); lde2
-    -%at('diff(y(t),t,1),t=0)+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s
+    (-%at('diff(y(t),t,1),t=0))+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s
 
 Dies besagt
 
diff --git a/src/doc/en/constructions/linear_algebra.rst b/src/doc/en/constructions/linear_algebra.rst
index 58747c663f1..e101cab2c01 100644
--- a/src/doc/en/constructions/linear_algebra.rst
+++ b/src/doc/en/constructions/linear_algebra.rst
@@ -417,12 +417,7 @@ Using maxima, you can easily solve linear equations:
     (a, b, c)
     sage: eqn = [a+b*c==1, b-a*c==0, a+b==5]
     sage: s = solve(eqn, a,b,c); s
-    [[a == (25*I*sqrt(79) + 25)/(6*I*sqrt(79) - 34),
-      b == (5*I*sqrt(79) + 5)/(I*sqrt(79) + 11),
-      c == 1/10*I*sqrt(79) + 1/10],
-     [a == (25*I*sqrt(79) - 25)/(6*I*sqrt(79) + 34),
-      b == (5*I*sqrt(79) - 5)/(I*sqrt(79) - 11),
-      c == -1/10*I*sqrt(79) + 1/10]]
+    [[a == 50/(I*sqrt(79) + 11), b == (5*I*sqrt(79) + 5)/(I*sqrt(79) + 11), c == 1/10*I*sqrt(79) + 1/10], [a == -50/(I*sqrt(79) - 11), b == (5*I*sqrt(79) - 5)/(I*sqrt(79) - 11), c == -1/10*I*sqrt(79) + 1/10]]
 
 You can even nicely typeset the solution in LaTeX:
 
diff --git a/src/doc/en/tutorial/tour_algebra.rst b/src/doc/en/tutorial/tour_algebra.rst
index a6b88e8301e..fa93dae888c 100644
--- a/src/doc/en/tutorial/tour_algebra.rst
+++ b/src/doc/en/tutorial/tour_algebra.rst
@@ -219,7 +219,7 @@ the notation :math:`x=x_{1}`, :math:`y=x_{2}`):
 
     sage: de1 = maxima("2*diff(x(t),t, 2) + 6*x(t) - 2*y(t)")
     sage: lde1 = de1.laplace("t","s"); lde1
-    2*(-%at('diff(x(t),t,1),t=0)+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)
+    2*((-%at('diff(x(t),t,1),t=0))+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)
 
 This is hard to read, but it says that
 
@@ -234,7 +234,7 @@ Laplace transform of the second equation:
 
     sage: de2 = maxima("diff(y(t),t, 2) + 2*y(t) - 2*x(t)")
     sage: lde2 = de2.laplace("t","s"); lde2
-    -%at('diff(y(t),t,1),t=0)+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s
+    (-%at('diff(y(t),t,1),t=0))+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s
 
 This says
 
diff --git a/src/doc/es/tutorial/tour_algebra.rst b/src/doc/es/tutorial/tour_algebra.rst
index 8f255eba1b1..e321d15f050 100644
--- a/src/doc/es/tutorial/tour_algebra.rst
+++ b/src/doc/es/tutorial/tour_algebra.rst
@@ -198,7 +198,7 @@ la notación :math:`x=x_{1}`, :math:`y=x_{2}`):
 
     sage: de1 = maxima("2*diff(x(t),t, 2) + 6*x(t) - 2*y(t)")
     sage: lde1 = de1.laplace("t","s"); lde1
-    2*(-%at('diff(x(t),t,1),t=0)+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)
+    2*((-%at('diff(x(t),t,1),t=0))+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)
 
 El resultado puede ser difícil de leer, pero significa que
 
@@ -213,7 +213,7 @@ Toma la transformada de Laplace de la segunda ecuación:
 
     sage: de2 = maxima("diff(y(t),t, 2) + 2*y(t) - 2*x(t)")
     sage: lde2 = de2.laplace("t","s"); lde2
-    -%at('diff(y(t),t,1),t=0)+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s
+    (-%at('diff(y(t),t,1),t=0))+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s
 
 Esto dice
 
diff --git a/src/doc/fr/tutorial/tour_algebra.rst b/src/doc/fr/tutorial/tour_algebra.rst
index e944975ed89..033af144450 100644
--- a/src/doc/fr/tutorial/tour_algebra.rst
+++ b/src/doc/fr/tutorial/tour_algebra.rst
@@ -183,7 +183,7 @@ Solution : Considérons la transformée de Laplace de la première équation
 
     sage: de1 = maxima("2*diff(x(t),t, 2) + 6*x(t) - 2*y(t)")
     sage: lde1 = de1.laplace("t","s"); lde1
-    2*(-%at('diff(x(t),t,1),t=0)+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)
+    2*((-%at('diff(x(t),t,1),t=0))+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)
 
 La réponse n'est pas très lisible, mais elle signifie que
 
@@ -198,7 +198,7 @@ la seconde équation :
 
     sage: de2 = maxima("diff(y(t),t, 2) + 2*y(t) - 2*x(t)")
     sage: lde2 = de2.laplace("t","s"); lde2
-    -%at('diff(y(t),t,1),t=0)+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s
+    (-%at('diff(y(t),t,1),t=0))+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s
 
 Ceci signifie
 
diff --git a/src/doc/ja/tutorial/tour_algebra.rst b/src/doc/ja/tutorial/tour_algebra.rst
index a82d16f31fa..04ac41e33b8 100644
--- a/src/doc/ja/tutorial/tour_algebra.rst
+++ b/src/doc/ja/tutorial/tour_algebra.rst
@@ -216,8 +216,7 @@ Sageを使って常微分方程式を研究することもできる． :math:`x'
 
     sage: de1 = maxima("2*diff(x(t),t, 2) + 6*x(t) - 2*y(t)")
     sage: lde1 = de1.laplace("t","s"); lde1
-    2*(-%at('diff(x(t),t,1),t=0)+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)
-
+    2*((-%at('diff(x(t),t,1),t=0))+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)
 
 この出力は読みにくいけれども，意味しているのは
 
@@ -231,7 +230,7 @@ Sageを使って常微分方程式を研究することもできる． :math:`x'
 
     sage: de2 = maxima("diff(y(t),t, 2) + 2*y(t) - 2*x(t)")
     sage: lde2 = de2.laplace("t","s"); lde2
-    -%at('diff(y(t),t,1),t=0)+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s
+    (-%at('diff(y(t),t,1),t=0))+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s
 
 意味するところは
 
diff --git a/src/doc/pt/tutorial/tour_algebra.rst b/src/doc/pt/tutorial/tour_algebra.rst
index f456b0f70e6..574661da790 100644
--- a/src/doc/pt/tutorial/tour_algebra.rst
+++ b/src/doc/pt/tutorial/tour_algebra.rst
@@ -207,7 +207,7 @@ equação (usando a notação :math:`x=x_{1}`, :math:`y=x_{2}`):
 
     sage: de1 = maxima("2*diff(x(t),t, 2) + 6*x(t) - 2*y(t)")
     sage: lde1 = de1.laplace("t","s"); lde1
-    2*(-%at('diff(x(t),t,1),t=0)+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)
+    2*((-%at('diff(x(t),t,1),t=0))+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)
 
 O resultado é um pouco difícil de ler, mas diz que
 
@@ -222,7 +222,7 @@ calcule a transformada de Laplace da segunda equação:
 
     sage: de2 = maxima("diff(y(t),t, 2) + 2*y(t) - 2*x(t)")
     sage: lde2 = de2.laplace("t","s"); lde2
-    -%at('diff(y(t),t,1),t=0)+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s
+    (-%at('diff(y(t),t,1),t=0))+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s
 
 O resultado significa que
 
diff --git a/src/doc/ru/tutorial/tour_algebra.rst b/src/doc/ru/tutorial/tour_algebra.rst
index fb2cac89968..30285e87113 100644
--- a/src/doc/ru/tutorial/tour_algebra.rst
+++ b/src/doc/ru/tutorial/tour_algebra.rst
@@ -200,7 +200,7 @@ Sage может использоваться для решения диффер
 
     sage: de1 = maxima("2*diff(x(t),t, 2) + 6*x(t) - 2*y(t)")
     sage: lde1 = de1.laplace("t","s"); lde1
-    2*(-%at('diff(x(t),t,1),t=0)+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)
+    2*((-%at('diff(x(t),t,1),t=0))+s^2*'laplace(x(t),t,s)-x(0)*s)-2*'laplace(y(t),t,s)+6*'laplace(x(t),t,s)
 
 Данный результат тяжело читаем, однако должен быть понят как
 
@@ -212,7 +212,7 @@ Sage может использоваться для решения диффер
 
     sage: de2 = maxima("diff(y(t),t, 2) + 2*y(t) - 2*x(t)")
     sage: lde2 = de2.laplace("t","s"); lde2
-    -%at('diff(y(t),t,1),t=0)+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s
+    (-%at('diff(y(t),t,1),t=0))+s^2*'laplace(y(t),t,s)+2*'laplace(y(t),t,s)-2*'laplace(x(t),t,s)-y(0)*s
 
 Результат:
 
diff --git a/src/sage/calculus/desolvers.py b/src/sage/calculus/desolvers.py
index e8e3277824a..1957bd6e75f 100644
--- a/src/sage/calculus/desolvers.py
+++ b/src/sage/calculus/desolvers.py
@@ -281,7 +281,7 @@ def desolve(de, dvar, ics=None, ivar=None, show_method=False, contrib_ode=False)
     Some more types of ODE's::
 
         sage: desolve(x*diff(y,x)^2-(1+x*y)*diff(y,x)+y==0,y,contrib_ode=True,show_method=True)
-        [[y(x) == _C*e^x, y(x) == _C + log(x)], 'factor']
+        [[y(x) == _C + log(x), y(x) == _C*e^x], 'factor']
 
     ::
 
diff --git a/src/sage/interfaces/maxima_abstract.py b/src/sage/interfaces/maxima_abstract.py
index 1392573fbad..e22f29c5d34 100644
--- a/src/sage/interfaces/maxima_abstract.py
+++ b/src/sage/interfaces/maxima_abstract.py
@@ -426,7 +426,7 @@ def version(self):
         EXAMPLES::
 
             sage: maxima.version()
-            '5.36.0.1'
+            '5.38.1'
         """
         return maxima_version()
 
@@ -2213,7 +2213,7 @@ def maxima_version():
 
         sage: from sage.interfaces.maxima_abstract import maxima_version
         sage: maxima_version()
-        '5.36.0.1'
+        '5.38.1'
     """
     return os.popen('maxima --version').read().split()[-1]
 
diff --git a/src/sage/interfaces/maxima_lib.py b/src/sage/interfaces/maxima_lib.py
index ed654b4655a..41ae822e4ff 100644
--- a/src/sage/interfaces/maxima_lib.py
+++ b/src/sage/interfaces/maxima_lib.py
@@ -1087,7 +1087,7 @@ def to_poly_solve(self,vars,options=""):
             sage: from sage.interfaces.maxima_lib import maxima_lib
             sage: sol = maxima_lib(sin(x) == 0).to_poly_solve(x)
             sage: sol.sage()
-            [[x == pi*z54]]
+            [[x == pi*z...]]
         """
         if options.find("use_grobner=true") != -1:
             cmd=EclObject([[max_to_poly_solve], self.ecl(), sr_to_max(vars),
diff --git a/src/sage/matrix/matrix1.pyx b/src/sage/matrix/matrix1.pyx
index f818644bf9f..86e13aace79 100644
--- a/src/sage/matrix/matrix1.pyx
+++ b/src/sage/matrix/matrix1.pyx
@@ -203,7 +203,7 @@ cdef class Matrix(matrix0.Matrix):
             sage: a = maxima(m); a
             matrix([0,1,2],[3,4,5],[6,7,8])
             sage: a.charpoly('x').expand()
-            -x^3+12*x^2+18*x
+            (-x^3)+12*x^2+18*x
             sage: m.charpoly()
             x^3 - 12*x^2 - 18*x
         """
diff --git a/src/sage/symbolic/constants.py b/src/sage/symbolic/constants.py
index 58ed87f289e..4b664926c8c 100644
--- a/src/sage/symbolic/constants.py
+++ b/src/sage/symbolic/constants.py
@@ -59,7 +59,7 @@
     sage: a = pi + e*4/5; a
     pi + 4/5*e
     sage: maxima(a)
-    %pi+4*%e/5
+    %pi+(4*%e)/5
     sage: RealField(15)(a)           # 15 *bits* of precision
     5.316
     sage: gp(a)
diff --git a/src/sage/symbolic/expression.pyx b/src/sage/symbolic/expression.pyx
index 5a46186d69e..044589f410f 100644
--- a/src/sage/symbolic/expression.pyx
+++ b/src/sage/symbolic/expression.pyx
@@ -11537,14 +11537,14 @@ cdef class Expression(CommutativeRingElement):
 
             sage: (n,k,j)=var('n,k,j')
             sage: sum(binomial(n,k)*binomial(k-1,j)*(-1)**(k-1-j),k,j+1,n)
-            -sum((-1)^(-j + k)*binomial(k - 1, j)*binomial(n, k), k, j + 1, n)
+            -(-1)^(-j)*sum((-1)^k*binomial(k - 1, j)*binomial(n, k), k, j + 1, n)
             sage: assume(j>-1)
             sage: sum(binomial(n,k)*binomial(k-1,j)*(-1)**(k-1-j),k,j+1,n)
             1
             sage: forget()
             sage: assume(n>=j)
             sage: sum(binomial(n,k)*binomial(k-1,j)*(-1)**(k-1-j),k,j+1,n)
-            -sum((-1)^(-j + k)*binomial(k - 1, j)*binomial(n, k), k, j + 1, n)
+            -(-1)^(-j)*sum((-1)^k*binomial(k - 1, j)*binomial(n, k), k, j + 1, n)
             sage: forget()
             sage: assume(j==-1)
             sage: sum(binomial(n,k)*binomial(k-1,j)*(-1)**(k-1-j),k,j+1,n)
@@ -11552,7 +11552,7 @@ cdef class Expression(CommutativeRingElement):
             sage: forget()
             sage: assume(j<-1)
             sage: sum(binomial(n,k)*binomial(k-1,j)*(-1)**(k-1-j),k,j+1,n)
-            -sum((-1)^(-j + k)*binomial(k - 1, j)*binomial(n, k), k, j + 1, n)
+            -(-1)^(-j)*sum((-1)^k*binomial(k - 1, j)*binomial(n, k), k, j + 1, n)
             sage: forget()
 
         Check that :trac:`16176` is fixed::
diff --git a/src/sage/tests/french_book/recequadiff.py b/src/sage/tests/french_book/recequadiff.py
index 687df8b0756..5ac6bfdfadf 100644
--- a/src/sage/tests/french_book/recequadiff.py
+++ b/src/sage/tests/french_book/recequadiff.py
@@ -195,14 +195,14 @@
 Sage example in ./recequadiff.tex, line 575::
 
   sage: Sol(x) = solve(sol, y)[0]; Sol(x)
-  log(y(x)) == (_C + x)*a + log(b*y(x) - a)
+  log(y(x)) == _C*a + a*x + log(b*y(x) - a)
 
 Sage example in ./recequadiff.tex, line 582::
 
   sage: Sol(x) = Sol(x).lhs()-Sol(x).rhs(); Sol(x)
-  -(_C + x)*a - log(b*y(x) - a) + log(y(x))
+  -_C*a - a*x - log(b*y(x) - a) + log(y(x))
   sage: Sol = Sol.simplify_log(); Sol(x)
-  -(_C + x)*a + log(y(x)/(b*y(x) - a))
+  -_C*a - a*x + log(y(x)/(b*y(x) - a))
   sage: solve(Sol, y)[0].simplify()
   y(x) == a*e^(_C*a + a*x)/(b*e^(_C*a + a*x) - 1)