As the user rburing advised in the thread https://ask.sagemath.org/question/44995/combine-plots-with-built-in-maxima-trajectory-in-sage-available/ I'm opening this one now.
When running the following code, one obtains a wrong output:
y=function('y')(x)
desolve(diff(y)==4*y/x+x*sqrt(y),y,ics=[1,1]).factor()
The output is 1/4*x^4*(log(x) - 2)^2
instead of 1/4*x^4*(log(x) + 2)^2
. Mathematica however outputs both (by running DSolve[{D[y[x], x] == 4*y[x]/x + x*Sqrt[y[x]], y[1] == 1}, y[x], x]
).