using SymPy

x = Sym("x")

@vars a b c
a,b,c = @syms a,b,c

u = symbols("u")
x = symbols("x", real=true)
y1, y2 = symbols("y1, y2", positive=true)
alpha = symbols("alpha", integer=true, positive=true)

solve(x^2 + 1)   # ±i are not real

solve(y1 + 1)    # -1 is not positive

PI, E, oo

asin(1), asin(Sym(1))

x, y = symbols("x,y")
ex = x^2 + 2x + 1
subs(ex, x, y)

subs(ex, x, 0)

x,y,z = symbols("x,y,z")
ex = x + y + z
subs(ex, (x,1), (y,pi))      

subs(ex, x=>1, y=>pi)

ex(x=>1, y=>pi)

ex(1, pi)

ex |> subs(x, 1)

ex |> replace(y, pi)

N(PI)  # floating-point value

evalf(PI)

N(PI, 30)

evalf(PI, 30)

x,y,z = symbols("x, y, z")

p = x^2 + 3x + 2
factor(p)

expand(prod([(x-i) for i in 1:5]))

factor(x^2 - 2)

q = x*y + x*y^2 + x^2*y + x

collect(q, x)

collect(q, y)

simplify(q)

expand((x + 1)*(x - 2) - (x - 1)*x)

factor(exp(2x) + 3exp(x) + 2)

r = 1/x + 1/x^2

together(r)

apart( (4x^3 + 21x^2 + 10x + 12) /  (x^4 + 5x^3 + 5x^2 + 4x))

top = (x-1)*(x-2)*(x-3)
bottom = (x-1)*(x-4)
top/bottom

r = expand(top) / expand(bottom)

cancel(r)

x,y = symbols("x,y", nonnegative=true)
a = symbols("a", real=true)
simplify(x^a * y^a - (x*y)^a)

x,y,a = symbols("x,y,a")
simplify(x^a * y^a - (x*y)^a)

powsimp(x^a * y^a - (x*y)^a, force=true)

theta = symbols("theta", real=true)
p = cos(theta)^2 + sin(theta)^2

simplify(p)

simplify(sin(2theta) - 2sin(theta)*cos(theta))

expand_trig(sin(2theta))

a,b,c,x = symbols("a, b, c, x") 
p = a*x^2 + b*x + c

coeff(p, x^2) # a 
coeff(p, x)   # b

subs(p, x, 0) # c

[[coeff(p, x^i) for i in N(degree(p)):-1:1]; subs(p,x,0)]

coeffs(p) # fails

q = Poly(p, x)
coeffs(q)

real_roots(x^2 - 2)

p = (x-3)^2*(x-2)*(x-1)*x*(x+1)*(x^2 + x + 1)
real_roots(p)

solve(p)

polyroots(p)

p = x^5 - x + 1
polyroots(p)

rts = solve(p)

[N(r) for r in rts]     # or map(N, rts)

nroots(p)

solve(cos(x) - sin(x))

solve(cos(x) - x)

a,b,c  = symbols("a,b,c", real=true)
p = a*x^2 + b*x + c
solve(p, x)

solve(p)

x, y = symbols("x,y", real=true)
exs = [2x+3y-6, 3x-4y-12]
d = solve(exs)

vars = [x,y]
vals = [d[string(var)] for var in vars]
[subs(ex, zip(vars, vals)...) for ex in exs]

a,b,c,h = symbols("a,b,c,h", real=true)
p = a*x^2 + b*x + c
fn = cos
exs = [fn(0*h)-subs(p,x,0), fn(h)-subs(p,x,h), fn(2h)-subs(p,x,2h)]
d = solve(exs, [a,b,c])

vars=[a,b,c]
vals = [d[string(var)] for var in vars]
quad_approx = subs(p, zip(vars, vals)...)

n = 3
x, c = symbols("x,c")
as = Sym["a$i" for i in 0:(n-1)]
bs = Sym["b$i" for i in 0:(n-1)]
p = sum([as[i+1]*x^i for i in 0:(n-1)])
q = sum([bs[i+1]*(x-c)^i for i in 0:(n-1)])
solve(p-q, bs)

solve(Eq(x, 1))

solve(x ⩵ 1)

x, y = symbols("x,y", real=true)
exs = [2x+3y ⩵ 6, 3x-4y ⩵ 12]    ## Using \Equal[tab]
d = solve(exs)

limit(sin(x)/x, x, 0)

limit((1+1/x)^x, x, oo)

a = symbols("a", positive=true)
ex = (sqrt(2a^3*x-x^4) - a*(a^2*x)^(1//3)) / (a - (a*x^3)^(1//4))

subs(ex, x, a)

subs(denom(ex), x, a), subs(numer(ex), x, a)

limit(ex, x, a)

limit(quad_approx, h, 0)

limit(sign(x), x, 0)

limit(sign(x), x, 0, dir="-"), limit(sign(x), x, 0, dir="+")

f(x) = sin(5x)/x
limit(f, 0)

f(x) = 1/x^(log(log(log(log(1/x)))) - 1)

hs = [10.0^(-i) for i in 6:16]
ys = [f(h) for h in hs]
[hs ys]

limit(f(x), x, 0, dir="+")

x, h = symbols("x,h")
f(x) = exp(x)*sin(x)
limit((f(x+h) - f(x)) / h, h, 0)

diff(f(x), x)

diff(x^x, x)

diff(exp(-x^2), x, 2)

diff(exp(-x^2), x, x, x)     # same as diff(..., x, 3)

f(x) = (12x^2 - 1) / (x^3)
diff(f(x), x) |> solve

f(x) = exp(x)*cos(x)
diff(f, 2)

x,y = symbols("x,y")
ex = x^2*cos(y)
Sym[diff(ex,v1, v2) for v1 in [x,y], v2 in [x,y]]

ex = Derivative(exp(x*y), x, y, 2)

doit(ex)

F, G = symbols("F, G", cls=SymFunction)

F(x)

diff(F(x))

x,y = symbols("x, y")
ex = y^4 - x^4 - y^2 + 2x^2

ex1 = subs(ex, y, F(x))

ex2 = diff(ex1, x)

ex3 = solve(ex2, diff(F(x)))[1]

ex4 = subs(ex3, F(x), y)

w, h, P = symbols("w, h, P", nonnegative=true)
r = w/2
A = w*h + 1//2 * (pi * r^2)
p = w + 2h + pi*r

h0 =  solve(P-p, h)[1]
A1 = subs(A, h, h0)

coeffs(Poly(A1, w))

coeff(collect(expand(A1), w), w^2)

subs(A1, w, 0)

b = solve(subs(P-p, h, 0), w)[1]

c = solve(diff(A1, w), w)[1]

atb = subs(A1, w, b)
atc = subs(A1, w, c)

atc - atb

simplify(atc - atb)

integrate(x^3, x)

n = symbols("n", real=true)
ex = integrate(x^n, x)

subs(ex, n, 3)

integrate(x^2, (x, 0, 1))

integrate(x^5*sin(x), x)

f(x) = (x^4 + x^2 * exp(x) - x^2 - 2x*exp(x) - 2x - exp(x)) * exp(x) / ( (x-1)^2 * (x+1)^2 * (exp(x) + 1) )
integrate(f(x), x)

x, y = symbols("x,y")
integrate(x*y, (y, 0, 1), (x, 0, 1))

integrate(x^2*y, (y, 0, sqrt(1 - x^2)), (x, -1, 1))

integ = Integral(sin(x^2), x)

doit(integ)

f(x) = exp(x) * cos(x)
integrate(f)

integrate(sin, 0, pi)

s1 = series(exp(sin(x)), x, 0, 4)

s2 = series(cos(exp(x)), x, 0, 6)

simplify(s1 * s2)

removeO(s1)

i, n = symbols("i, n")
summation(i^2, (i, 1, n))

sn = Sum(1/i^2, (i, 1, n))
doit(sn)

limit(doit(sn), n, oo)

x,y = symbols("x,y")
v = [1,2,x]
w = [1,y,3]

dot(v,w)

cross(v,w)

ex = x^2*y - x*y^2
Sym[diff(ex,var) for var in [x,y]]

Sym[diff(ex, v1, v2) for v1 in [x,y], v2 in [x,y]]

hessian(ex)

x,y = symbols("x,y")
M = [1 x; x 1]

diagm(ones(Sym, 5))
M^2
det(M)

rref(M)

M[:rref]()

M = [one(Sym) 2 3 0 0; 4 10 0 0 1]
vs = nullspace(M)

[M*vs[i] for i in 1:3]

M = [1 x; x 1]
P, D = diagonalize(M)  # M = PDP^-1
D

F = symbols("F", cls=SymFunction)

diffeq = Eq(diff(F(x), x, 2) - 2*diff(F(x)) + F(x), sin(x))

ex = dsolve(diffeq, F(x))

ex1 = rhs(ex)

solve(subs(ex1, x, 0), Sym("C1"))

ex2 = subs(ex1, Sym("C1"), -1//2)

solve( subs(diff(ex2, x), x, 0) - 1, Sym("C2") )

ex3 = subs(ex2, Sym("C2"), 3//2)

t, m,k,alpha = symbols("t,m,k,alpha")
v = symbols("v", cls=SymFunction)
ex = Eq( (m/k)*diff(v(t),t), alpha^2 - v(t)^2 )

ex[:classify_ode]()

dsolve(ex, v(t))