Lorenz_RK4_BT.jl

function lorenz(h,N,σ,ρ,β,x0,y0,z0,a,b,c)
    # stage
    s = length(b)
    # Number of Equations
    M = 3
    #--------------------
    # Equations
    #--------------------
    f₁(x,y,z) = -σ * x + σ * y
    f₂(x,y,z) = -x * z + ρ * x - y
    f₃(x,y,z) =  x * y - β * z

    f(x) = [f₁(x...) f₂(x...) f₃(x...)]
    #--------------------
    # Initialization
    #--------------------
    x = zeros(N,M)
    k = zeros(N,s,M)
    #----------------------
    # Initial Value
    #----------------------
    x[1,1] = x0
    x[1,2] = y0
    x[1,3] = z0
    #----------------------
    # Numerical Calculation
    #----------------------
    @time for n = 1:N-1
        # k
        for i = 1:s
            k[n,i,:] = f( x[n:n,:] + h * a[i:i, 1:i-1] * k[n, 1:i-1, :] )
        end
        # x
        x[n+1:n+1,:] = x[n:n,:] + h * b * k[n,:,:]
    end
    return x
end

#----------------------
# Butcher Tableau (RK4)
#----------------------
a = [ 0   0   0   0
      1/2 0   0   0
      0   1/2 0   0
      0   0   1   0 ]

b = [0 1/2 1/2 1]

c = [1/6 1/3 1/3 1/6]

#----------------------
# Conditions
#----------------------
# step
h = 0.001
# Number of Points
N  = 30000
# limit
tlm= [0,(N-1)*h]
# t 
t  = tlm[1]:h:tlm[2]
#----------------------
# Parameter
#----------------------
σ = 10.
ρ = 28.
β = 8.0 / 3.0
#----------------------
# Initial Value
#----------------------
x0 = 0.
y0 = 1.
z0 = 1.05
#----------------------
# Numerical Calculation
#----------------------
x = lorenz(h,N,σ,ρ,β,x0,y0,z0,a,b,c)
#----------------------
# Plot
#----------------------
using Plots
gr()

plot(x[:,1],x[:,2],x[:,3])
savefig("Lorenz_RK4_BT.png")