using DifferentialEquations
using Plots

# Simulate Lorenz 63 system and investigate sensitivity to initial conditions

function tendency!(du, u, p, t)
    x,y,z = u
    σ,ρ,β = p
    
    du[1] = dx = σ*(y-x)
    du[2] = dy = x*(ρ-z) - y
    du[3] = dz = x*y - β*z
end

p = [10.0, 28.0, 8/3] # parameters of the Lorentz 63 system
tspan = (0.0, 1000.0)

u0 = [1.0, 0.0, 0.0]
δu0 = [0.0, 0.001, 0.0] # small deviation in initial condition

prob = ODEProblem(tendency!, u0, tspan, p)
sol1 = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8) # control simulation


prob = ODEProblem(tendency!, u0 + δu0, tspan, p)
sol2 = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8) # perturbed initial condition

plot(sol1, idxs = (1, 2, 3)) # 3d plot

# plot(sol1, idxs = (1, 2))

plot([sol1[1, :], sol2[1, :]], [sol1[2, :], sol2[2, :]])