Runge-Kutta

pkg"add StaticArrays"
#  starting node
using StaticArrays, BenchmarkTools

const Vec3 = SVector{3, Float64}

function f(t, X, p)
    Vec3(
        p.a - p.b * X[1] + X[3] - p.d * X[1], 
        p.b * X[1] * X[3] - p.s * X[2], 
        p.m * X[2]-(p.p1 + p.p2) * X[3] - p.b * X[1] * X[3]
    )
end

function main()
    x0 = 0
    #ending node
    T = 1000
    params = (
        a = 0.056, b = 0.001, d = 0.0024, s = 0.03, m = 0.013, p1 = 1, p2 = 4
    )
    IC = Vec3(400, 0, 10)
    h = 0.01
    N = round(Int, abs(T-x0) / h)
    y = fill(zero(IC), N + 1)
    y[1] = IC
    t = x0:h:T
    # Body of Algorithm
    for i = 1:N
        K1 = h .* f(t[i], y[i], params)
        K2 = h .* f(t[i] + 0.5 * h, y[i] .+ K1 ./ 2.0, params)
        K3 = h .* f(t[i] + 0.5 * h, y[i] .+ K2 ./ 2.0, params)
        K4 = h .* f(t[i] + h, y[i] .+ K3, params)
        kk = (K1 .+ K2 .* 2.0 .+ K3 .* 2.0 .+ K4)
        y[i + 1] = y[i] .+ kk ./ 6.0
    end
end
@btime main()