function myfft(x)
    N=length(x)
    if N == 1
        return x
    end
    if N % 2 == 1
        println("N = ",N," was not even!")
        throw(ArgumentError)
    end
    ye=myfft(x[1:2:N-1])
    yo=myfft(x[2:2:N])
    M = N ÷ 2
    y = Array{Complex}(undef,N)
    for l = 1:M
        z0 = ye[l]
        z1 = exp(-1im*2*pi*(l-1)/N)*yo[l]
        y[l] = z0 + z1
        y[l+M] = z0 - z1
    end
    return y
end

function myifftwork(x)
    N=length(x)
    if N == 1
        return x
    end
    if N % 2 == 1
        println("N = ",N," was not even!")
        throw(ArgumentError)
    end
    ye=myifftwork(x[1:2:N-1])
    yo=myifftwork(x[2:2:N])
    M = N ÷ 2
    y = Array{Complex}(undef,N)
    for l = 1:M
        z0 = ye[l]
        z1 = exp(1im*2*pi*(l-1)/N)*yo[l]
        y[l] = z0 + z1
        y[l+M] = z0 - z1
    end
    return y
end

function myifft(x)
    return myifftwork(x)/length(x)
end

function f(x)
    if x>1 && x<2
        return 1
    else
        return 0
    end
end

function g(t,b)
    return -F.*b
end
function euler(t,b)
    return b+dT*g(t,b)
end
function rk4(t,b)
    k1=g(t,b)
    k2=g(t+dT/2,b+dT*k1/2)
    k3=g(t+dT/2,b+dT*k2/2)
    k4=g(t+dT,b+dT*k3)
    return b+(1/6)*dT*(k1+2*k2+2*k3+k4)
end

using Plots

L=2*pi
N=64
h=L/N
X=0:h:L-h
Y0=f.(X)
B0=myifft(Y0)
K=N÷2-1
D=zeros(Complex,N)
D[1:K+1]=0:K
D[N-K+1:N]=-K:-1
F=D.*D*(2*pi/L)^2

T=0.5
Tsteps=1024
dT=T/Tsteps

Bn=B0
for i=0:Tsteps-1
    global Bn
    Tn=i*dT
    Bn=rk4(Tn,Bn)
    if i%16==0
        Yn=real.(myfft(Bn))
        display(plot(X,Y0))
        display(plot!(X,Yn))
    end
end
Yn=real.(myfft(Bn))
display(plot(X,Y0))
display(plot!(X,Yn))