lorenz

function lorenz(h,n,a,b,c,x0,y0,z0)
    #--------------------
    # Equation
    #--------------------
    f₁(x,y,z) = -a * x + a * y
    f₂(x,y,z) = -x * z + b * x - y
    f₃(x,y,z) =  x * y - c * z
    #--------------------
    # Initialize
    #--------------------
    x = zeros(n)
    y = zeros(n)
    z = zeros(n)
    #----------------------
    # Initial Value
    #----------------------
    x[1] = x0
    y[1] = y0
    z[1] = z0
    #----------------------
    # Numerical Calculation
    #----------------------
    @time for i in 1:n-1
        kx₁ = f₁( x[i], y[i], z[i] )
        ky₁ = f₂( x[i], y[i], z[i] )
        kz₁ = f₃( x[i], y[i], z[i] )

        kx₂ = f₁( x[i] + h / 2. * kx₁, y[i] + h / 2. * ky₁, z[i] + h / 2. * kz₁)
        ky₂ = f₂( x[i] + h / 2. * kx₁, y[i] + h / 2. * ky₁, z[i] + h / 2. * kz₁)
        kz₂ = f₃( x[i] + h / 2. * kx₁, y[i] + h / 2. * ky₁, z[i] + h / 2. * kz₁)

        kx₃ = f₁( x[i] + h / 2. * kx₂, y[i] + h / 2. * ky₂, z[i] + h / 2. * kz₂)
        ky₃ = f₂( x[i] + h / 2. * kx₂, y[i] + h / 2. * ky₂, z[i] + h / 2. * kz₂)
        kz₃ = f₃( x[i] + h / 2. * kx₂, y[i] + h / 2. * ky₂, z[i] + h / 2. * kz₂)

        kx₄ = f₁( x[i] + h * kx₃, y[i] + h *  ky₃,　z[i] + h * kz₃)
        ky₄ = f₂( x[i] + h * kx₃, y[i] + h *  ky₃,　z[i] + h * kz₃)
        kz₄ = f₃( x[i] + h * kx₃, y[i] + h *  ky₃,　z[i] + h * kz₃)

        x[i+1] = x[i] + h / 6. * (kx₁ + 2. * kx₂ + 2. * kx₃ + kx₄)
        y[i+1] = y[i] + h / 6. * (ky₁ + 2. * ky₂ + 2. * ky₃ + ky₄)
        z[i+1] = z[i] + h / 6. * (kz₁ + 2. * kz₂ + 2. * kz₃ + kz₄)
    end
    return x,y,z
end
#----------------------
# Conditions
#----------------------
# Time Step [s]
dt = 0.001
# Number
n  = 30000
# limit [s]
tlm= [0,(n-1)*dt]
# t [s]
t  = tlm[1]:dt:tlm[2]
#----------------------
# Parameter
#----------------------
a = 10.
b = 28.
c = 8.0 / 3.0
#----------------------
# Initial Value
#----------------------
x0 = 0.
y0 = 1.
z0 = 1.05
#----------------------
# Numerical Calculation
#----------------------
x,y,z = lorenz(dt,n,a,b,c,x0,y0,z0)
#----------------------
# Plot
#----------------------
using Plots
gr()

plot(x,y,z)
savefig("Lorenz.png")