using PyPlot
x = linspace(0,1)
uref(x) = -x.^2/2 + x/2
plot(x, uref(x))

using Sundials

v0 = 0.0
function ode_rhs(t, uv, uvdot) 
    uvdot[1] = uv[2]
    uvdot[2] = -1.0
end
res = Sundials.cvode(ode_rhs, [0.0; v0], x)
plot(x, uref(x))
plot(x, res[:,1], "--")

plot(x, uref(x))
function ufinal(v0)
    res = Sundials.cvode(ode_rhs, [0.0; v0], x)
    plot(x, res[:,1], "--")
    return res[end,1]
end

function bisection_solver(f, xlo, xhi, tol=1e-3)
    flo = f(xlo)
    fhi = f(xhi)
    print("Initial bracket: [$xlo, $xhi]\n")
    while abs(xlo-xhi) > tol
        xmid = (xlo+xhi)/2
        fmid = f(xmid)
        if flo*fmid > 0.0
            xlo, flo = xmid, fmid
        elseif flo*fmid < 0.0
            xhi, fhi = xmid, fmid
        end
        print("  refine: [$xlo, $xhi]\n")
    end
    (xlo, xhi)
end

bisection_solver(ufinal, 0.0, 2.0)

NN = length(x)
N = NN-2
Tdense(N) = 2*eye(N) - diagm(ones(N-1),1) - diagm(ones(N-1),-1)
T = Tdense(N)
T[1:5,1:5]

h = 1.0/(N+1)
ufd = zeros(NN)
f = ones(NN)
ufd[2:N+1] = T\(h^2*f[2:N+1])
plot(x,uref(x))
plot(x,ufd, "--")

Tsparse(N) = 2*speye(N) - spdiagm(ones(N-1),1,N,N) - spdiagm(ones(N-1),-1,N,N)
T = Tsparse(N)
T[1:5,1:5]

ufd[2:N+1] = T\(h^2*f[2:N+1])
plot(x, uref(x))
plot(x, ufd, "--")

plot(x, uref(x)-ufd)

uref(x) = sin(pi*x)
fref(x) = pi^2 * sin(pi*x)

mmax = 8
hs = zeros(mmax)
errnorm = zeros(mmax)
for m = 3:mmax
    h = 2.0^-m
    N = 2^m-1
    x = linspace(0,1,N+2)
    ufd = Tsparse(N)\(h^2*fref(x[2:end-1]))
    hs[m] = h
    errnorm[m] = norm(uref(x[2:end-1])-ufd, Inf)
end
loglog(hs[3:mmax], errnorm[3:mmax])

errnorm[4:mmax]./errnorm[3:mmax-1]

N = 20
h = 1.0/(N+1)
x = linspace(0,1,N+2)

T = Tsparse(N)
u = zeros(N)
f = ones(N)

kappa = 1
dt = 1e-2
usteady(x) = (-x.^2+x)/2

plot(x[2:end-1], u)
for j = 1:10
    u += (-(kappa/h^2)*(T*u) + f) * dt
    plot(x[2:end-1], u)
end

(D,V) = eig(Tdense(100))
lambdas = D
print(lambdas[1], " ", lambdas[end], "\n")

u = zeros(N)
plot(x[2:end-1], u)
A = speye(N) + dt*(kappa/h^2)*T
for j = 1:50
    u = A\(u+dt*f)
    plot(x[2:end-1], u)
end
plot(x, usteady(x), "k--")

u = zeros(N)
A0 = speye(N) + dt*(kappa/h^2/2)*T
A1 = speye(N) - dt*(kappa/h^2/2)*T
for j = 1:50
    u = A0\(A1*u+dt*f)
    plot(x[2:end-1], u)
end
plot(x, usteady(x), "k--")

N = 100
x = linspace(0,1,N+2)
u0 = exp(-100*(x-0.5).^2)
u0 -= u0[1]
plot(x,u0)

xi = 0.4
u0 = exp(-100*(x-0.5).^2)
u0 -= u0[1]
u1 = copy(u0)
unext = copy(u1)
T = Tsparse(N)
for j = 1:100
    unext[2:end-1] = 2*u1[2:end-1]-u0[2:end-1]-xi^2*(T*u1[2:end-1])
    if mod(j,10) == 0
        plot(x, unext)
    end
    u0, u1, unext = u1, unext, u0
end