using Plots, ComplexPhasePortrait, ApproxFun, SingularIntegralEquations, DifferentialEquations
gr();

f = Fun(x -> cos(x^2), Chebyshev())   # f is expanded in Chebyshev coefficients
n = ncoefficients(f)   # This is the number of coefficients
D = zeros(n-1,n)
for k=1:n-1
    D[k,k+1] = k
end
D

fp = Fun(Ultraspherical(1), D*f.coefficients)

fp(0.1)

f'(0.1)

-2*0.1*sin(0.1^2)

D = Derivative()
(D*f)(0.1)

D : Chebyshev() → Ultraspherical(1) 

S₀ = I : Chebyshev() → Ultraspherical(1)

f = Fun(exp, Chebyshev())
g = S₀*f

g(0.1) - exp(0.1)

Multiplication(Fun(), Ultraspherical(2))'

B  = Evaluation(0.0) : Chebyshev()
D  = Derivative() : Chebyshev() → Ultraspherical(1)
S₀ = I : Chebyshev() → Ultraspherical(1)
L = [B; 
     D - S₀]

u = L \ [1; 0]
plot(u)

u(0.1) - exp(0.1)

D₀ = Derivative() : Chebyshev() → Ultraspherical(1)
D₁ = Derivative() : Ultraspherical(1) → Ultraspherical(2)
D₁*D₀  # 2 zeros not printed in (1,1) and (1,2) entry

S₀ = I : Chebyshev() → Ultraspherical(1)
S₁ = I : Ultraspherical(1) → Ultraspherical(2)
S₁*S₀

S₁*D₀  # or could have been D₁*S₀

B₋₁ = Evaluation(-1) : Chebyshev()
B₁ = Evaluation(1) : Chebyshev()
# u(-1) 
# u(1)
# u'' + u' +u

L = [B₋₁; 
     B₁; 
     D₁*D₀ + S₁*D₀ + S₁*S₀]

u = L \ [1.0,0.0,0.0]
plot(u)

x = Fun()
Jᵗ = Multiplication(x) : Chebyshev() → Chebyshev()  # transpose of the Jacobi operator

L = [B₋₁;   # u(-1)
     B₁ ;   # u(1)
     D₁*D₀ - S₁*S₀*Jᵗ]   # u'' - x*u

u = L \ [1.0;0.0;0.0]
plot(u; legend=false)

ε = 1E-6
L = [B₋₁; 
     B₁ ; 
     ε*D₁*D₀ - S₁*S₀*Jᵗ]

u = L \ [1.0;0.0;0.0]
plot(u; legend=false)

ε = 1E-10
L = [B₋₁; 
     B₁ ; 
     ε*D₁*D₀ - S₁*S₀*Jᵗ]


@time u = L \ [1.0;0.0;0.0]
@show ncoefficients(u);

M = -I + Jᵗ + (Jᵗ)^2  # represents -1+x+x^2

ε = 1E-6
L = [B₋₁; 
     B₁ ; 
     ε*D₁*D₀ - S₁*S₀*M]

@time u = L \ [1.0;0.0;0.0]

@show ε*u''(0.1) - (-1+0.1+0.1^2)*u(0.1)
plot(u)

p = Fun(exp, Chebyshev()) # polynomial approximation to exp(x)
M = Multiplication(p) : Chebyshev() # constructed using Clenshaw:

ApproxFun.bandwidths(M) # still banded

ε = 1E-6
L = [B₋₁; 
     B₁ ; 
     ε*D₁*D₀ + S₁*S₀*M]

@time u = L \ [1.0;0.0;0.0]

@show ε*u''(0.1) + exp(0.1)*u(0.1)
plot(u)