I forgot to mention SymPy. It is a symbolic toolbox, which cannot nearly match Maxima, but it aims to. Sage uses Maxima from Python, which might be a better idea.

Here is a session with SymPy in the terminal of Python.

>>> from sympy import *
>>> x=Symbol('x')
>>> expand((1+x)**10)
x**10 + 10*x**9 + 45*x**8 + 120*x**7 + 210*x**6 + 
   252*x**5 + 210*x**4 + 120*x**3 + 45*x**2 + 10*x + 1
>>> factor(diff(expand(1+x)**10,x))
10*(x + 1)**9

Many things like partial fractions, differentiation, integrals etc. work. But there is a lot of work to be done to get a useful production system.

The package seems to have serious bugs too.

>>> print limit((1+1/x)**x,x,oo)
E
>>> print limit(1+(1+1/x)**x,x,oo)
2

It does also have plotting features. The default backend is MatPlotLib.

from sympy import *
from sympy.plotting import *
u,v = symbols('u v')

plot3d_parametric_surface(u**3, v, u**2+v**2, (u, -5, 5), (v, -5, 5))

will produce the following plot.

As you see, the plots scale to a square plot automatically. This is something, Euler does not do, since the idea of Euler plots in 3D is a geometric view. But you can do it.

>plot3d("x^3","y","x^2+y^2",r=5,>spectral,scale=[1,25,5], ...
   angle=30°,height=30°,n=100,grid=20):

Update

I asked the people in the Google group „sympy“ about the limit bug. It is due to the heuristic

\(1^\infty = 1\)

which is implied in cases, where it does not hold. While the first limit is known, the known facts are not detected in the second limit, where 1+ is in front of the known fact.  They say, that internally it is difficult to change the behavior of operators inside limit. But that is exactly, what is necessary. (By the way, it is the same error, students often make.)

There is a more advanced limit method „gruntz“, which uses Taylor series. This method is the same as Maxima’s tlimit, which is now incorporated in Maxima’s limit. It can solve the limit and similar cases.

>>> from sympy.series import gruntz
>>> gruntz(1+(1+1/x)**x,x,oo)
1 + E
>>> gruntz((1+1/x)**x-E,x,oo)
0
>>> gruntz(((1+1/x)**x-E)*x,x,oo)
-E/2
>>> gruntz((((1+1/x)**x-E)*x+E/2)*x,x,oo)
11*E/24