SageMath

x,a,b,c = var('x,a,b,c')
log(sqrt(a)).simplify_log() # returns log(a)/2
sin(a+b).simplify_trig() # returns cos(a)*sin(b) + sin(a)*cos(b)
cos(a+b).simplify_trig() # returns cos(a)*cos(b) - sin(a)*sin(b)

limit((xˆ2+1)/(2+x+3*xˆ2), x=infinity) # returns 1/3
limit(sin(x)/x, x=0) # returns 1

diff(acos(x),x) # returns -1/sqrt(1 - xˆ2)
f = exp(x)*log(x)
f.diff(x,3) # returns e^x*log(x) + 3*e^x/x - 3*e^x/x^2 + 2*e^x/x^3

solve(a*x^2 + b*x + c, x) # returns [x == (-sqrt(b^2 - 4*a*c) - b)/(2*a),
 # x == (sqrt(b^2 - 4*a*c) - b)/(2*a)]

t = var('t') # define a variable t
x = function('x',t) # define x to be a function of that variable
DE = lambda y: diff(y,t) + y - 1
desolve(DE(x(t)), [x,t]) # returns '%e^-t*(%e^t+%c)'

A = Matrix([[1,2,3],[3,2,1],[1,1,1]])
y = vector([0,-4,-1])
A.solve_right(y) # returns (-2, 1, 0)
A.eigenvalues() # returns [5, 0, -1]

B = Matrix([[1,2,3],[3,2,1],[1,2,1]])
B.inverse() # returns [ 0 1/2 -1/2]
 # [-1/4 -1/4 1]
 # [ 1/2 0 -1/2]

# Call numpy for the Moore-Penrose pseudo-inverse,
# since Sage does not support that yet.

import numpy
C = Matrix([[1 , 1], [2 , 2]])
matrix(numpy.linalg.pinv(C.numpy())) # returns [0.1 0.2]
 # [0.1 0.2]

prime_pi(1000000) # returns 78498, the number of primes less than one million

E = EllipticCurve('389a') # construct an elliptic curve from its Cremona label
P, Q = E.gens()
7*P + Q # returns (2869/676 : -171989/17576 : 1)