#  rossler.jl -- integration of a chaotic non-linear system

using Printf, LinearAlgebra, Plots

a=0.2; b=0.2; c=5.7
logYs=zeros(10); loghs=zeros(10)

function f(y)
    r=[-y[2]-y[3],
       y[1]+a*y[2],
       b+y[3]*(y[1]-c)]
    return r
end

function euler(t,y,h)
    return y+h*f(y)
end

function rk2(t,y,h)
    k1=f(y)
    k2=f(y+h*2/3*k1)
    return y+h*(1/4*k1+3/4*k2)
end

function solve(t0,y0,h,n)
    yn=copy(y0)
    for j=1:n
        tn=t0+(j-1)*h
        yn=rk2(tn,yn,h)
    end
    return yn
end

function main()
    t0=0; T=1
    y0=[1.0,1,1]
    yold=zeros(3)
    @printf("%5s %7s %14s    %14s %14s %14s\n",
        "p","n","h","y1","y2","y3")
    for p=6:16
        n=2^p
        h=T/n
        yn=solve(t0,y0,h,n)
        @printf("%5d %7d %14.7e    %14.7e %14.7e %14.7e\n",
            p,n,h,yn[1],yn[2],yn[3])
        if p>6
            t=log(norm(yold-yn))
            global logYs[p-6]=t
            global loghs[p-6]=log(2*h)
        end
        yold=copy(yn)
    end
end

main()