diff --git a/examples/henon_heiles.jl b/examples/henon_heiles.jl
new file mode 100644
index 0000000..6dfeb0a
--- /dev/null
+++ b/examples/henon_heiles.jl
@@ -0,0 +1,21 @@
+using Distributed
+@everywhere using DifferentialEquations: ODEProblem, solve, ImplicitMidpoint
+
+struct HenonHeiles end
+
+@everywhere function deterministic_system(model::HenonHeiles)
+    function tendency!(dudt,u,p,t)
+        dudt[1] = u[3]
+        dudt[2] = u[4]
+        dudt[3] = -u[1]-2.0*u[1]*u[2]
+        dudt[4] = -u[2]-(u[1]^2-u[2]^2)
+    end
+    return tendency!
+end
+
+@everywhere function evolve_deterministic_system(model::HenonHeiles, u0, tspan, dt, alg_kwargs)
+    tendency! = deterministic_system(model)
+    prob = ODEProblem(tendency!, u0, tspan, nothing)
+    sol = solve(prob, ImplicitMidpoint(), dt = dt, saveat = tspan[1]:dt:tspan[2]; alg_kwargs...);
+    return sol.u
+end
diff --git a/examples/henon_heiles_integration_test.jl b/examples/henon_heiles_integration_test.jl
new file mode 100644
index 0000000..ee635cc
--- /dev/null
+++ b/examples/henon_heiles_integration_test.jl
@@ -0,0 +1,124 @@
+using Revise
+using RareEvents
+using Statistics
+using DelimitedFiles
+using LinearAlgebra
+using SpecialFunctions
+using Plots
+
+examples_dir = joinpath(pkgdir(RareEvents), "examples")
+include(joinpath(examples_dir, "henon_heiles.jl"))
+n_processes = Sys.CPU_THREADS
+addprocs(n_processes)
+
+# Set up model and wrapper function that evolves a single trajectory
+dt = 1.0
+alg_kwargs = ();
+FT = Float64
+model = HenonHeiles()
+@everywhere evolve_wrapper(tspan) = (u0) -> evolve_deterministic_system(model, u0, tspan, dt, alg_kwargs)
+
+
+# Set up rare event algorithm
+τ = 1.0
+function set_ic_via_y(E, y; py = 0.0, x = 0.0)
+    twiceV = (x^2 + y^2 + 2 * x^2 * y - 2 / 3 * y^3)
+    px = sqrt(2.0 * E - py^2 - twiceV)
+    return [x,y, px, py]
+end;
+E = 0.125;
+yvals = -0.35:0.001:0.35;
+nensemble = length(yvals)
+u0 = [set_ic_via_y(E, y) for y in yvals]
+k = 0.0
+ϵ = 0.000
+metric(y) = y[2]
+t_end = 1000.0
+t_start = 0.0
+tspan = (t_start, t_end)
+# averaging window
+T = 50
+steps_per_window = Int64.(round(T/dt))
+
+
+# Run GKLT -  bootstrap iterations - and compute necessary statistics per iteration
+iters = 5
+em = zeros(nensemble*iters)
+em = reshape(em, (nensemble,iters))
+r = similar(em)
+r_paper = similar(em)
+σr = similar(em)
+bin = 0.03
+a_range = Array(0.0:bin:0.6)
+counts = zeros(length(a_range))
+algorithm_cost = nensemble*(t_end-t_start)*iters*dt
+for iter in 1:iters
+    println(iter)
+    sim = RareEventSampler{Float64}(dt, τ, tspan,u0, evolve_wrapper, nensemble, metric, k, ϵ, iter)
+    
+    run!(sim);
+    a_m = zeros(nensemble)
+    lr_vector = zeros(nensemble)
+    λ = sum(log.(sim.weight_norm))/(t_end-t_start)
+    
+    for i in 1:nensemble
+        trajectory = [metric(sim.ensemble[i][j]) for j in 1:length(sim.ensemble[i])]
+        lr_vector[i] = likelihood_ratio(sim, trajectory, λ, t_end-t_start)
+        # This just counts the number of (treated as independent) block means per bin, over all trajectories, over all iterations
+        counts .+=  count_events_per_bin(block_mean(trajectory, steps_per_window), lr_vector[i],a_range)
+        # This is for estimating return times via the Ragone paper method; via `return_curve`.
+        a_m[i] = maximum(block_mean(trajectory, steps_per_window))
+
+    end
+    em[:,iter], r[:,iter], r_paper[:,iter], σr[:, iter]= return_curve(a_m, t_end-t_start, lr_vector);
+
+end
+rmprocs([p for p in procs() if p != myid()])
+
+
+# Evolve the solution directly; compute statistics
+tspan_direct = (0.0, 3.e6)
+direct_solution = evolve_deterministic_system(model, u0[end], tspan_direct, dt, (maxiters = 1e8,));
+direct_cost = (tspan_direct[2]-tspan_direct[1])*dt
+u = [metric(direct_solution[k]) for k in 1:length(direct_solution)]
+A = block_mean(u, steps_per_window)
+m = 100
+em_direct, r_direct, r_paper_direct, σr_direct = return_curve(A, FT(T),  FT.(ones(length(A))))
+
+# Fit GEV
+m = 100 
+blocks = block_maxima(A, m)
+params, nll = fit_gev(blocks, [std(blocks), mean(blocks), 0.1])
+cdf_A= (gev_cdf.(a_range, params[1], params[2], params[3])).^(1.0/m)
+pdf_A= (cdf_A[2:end] .- cdf_A[1:end-1])/bin
+
+
+# Make plots
+plot1 = plot()
+mean_a = zeros(length(a_range)-1)
+mean_rtn = zeros(length(a_range)-1)
+std_rtn = zeros(length(a_range)-1)
+for i in 1:(length(a_range)-1)
+    mask = (em[:] .< a_range[i+1]) .& ( em[:] .>= a_range[i])
+    if sum(mask) >0
+        mean_a[i] = mean(em[:][mask])
+        mean_rtn[i] = log10.(mean(r_paper[:][mask]))
+        std_rtn[i] = std(r_paper[:][mask])/mean(r_paper[:][mask])./log(10.0)
+    end
+end
+nonzero = (mean_a .!= 0.0) .& (isnan.(std_rtn) .== 0)
+final_a = mean_a[nonzero]
+final_r = mean_rtn[nonzero]
+final_σr = std_rtn[nonzero]
+plot!(plot1,final_a, final_r, ribbon = final_σr, label = "GKLT Algorithm; using maxes")
+denominator = 1.0 .- cumsum(counts)/sum(counts)
+safe_denominator = denominator .> 10.0*eps(FT)
+plot!(plot1,a_range[safe_denominator], log10.(50.0 ./denominator[safe_denominator]), label = "GKLT Algorithm; using counts")
+plot!(plot1,em_direct, log10.(r_paper_direct), ribbon = σr_direct ./ r_paper_direct ./ log(10.0), label = "Direct Integration")
+plot!(plot1, a_range, log10.(50.0 ./ (1.0 .- cdf_A)), label = "GEV")
+plot!(plot1, a_range, log10.(50.0 ./ (1.0 .- 0.5 .*(1 .+erf.(a_range ./2 .^0.5 ./std(A))))), label = "Gaussian")
+plot!(plot1,xlabel = "Event magnitude")
+plot!(plot1,ylabel = "Log10(Return Time)")
+plot!(plot1, legend = :topleft)
+plot!(plot1, xrange = [0.0,1.0])
+savefig(plot1, "return_curve_ou.png")
diff --git a/examples/standard_map.jl b/examples/standard_map.jl
new file mode 100644
index 0000000..137db50
--- /dev/null
+++ b/examples/standard_map.jl
@@ -0,0 +1,23 @@
+struct StandardMap{FT<: AbstractFloat}
+    K::FT
+end
+
+function mapping_system(model::StandardMap)
+    function update!(u)
+        u[4] = u[4] + model.K*cos(u[1])*u[3] # δI
+        u[3] = u[3] + u[4] # δθ
+        u[2] = mod(u[2] + model.K*sin(u[1]),2π) # I
+        u[1] = mod(u[1] + u[2],2π) # θ
+    end
+    return update!
+end
+
+function evolve_mapping_system(model::StandardMap, u0, tspan)
+    update! = mapping_system(model)
+    nsteps = Int(tspan[2]-tspan[1])
+    solution = map(0:nsteps) do (i)
+        update!(u0)
+        copy(u0)
+    end
+    return solution
+end
diff --git a/examples/standard_map_integration_test.jl b/examples/standard_map_integration_test.jl
new file mode 100644
index 0000000..99598ec
--- /dev/null
+++ b/examples/standard_map_integration_test.jl
@@ -0,0 +1,55 @@
+using Revise
+using RareEvents
+using Statistics
+using DelimitedFiles
+using Plots
+
+examples_dir = joinpath(pkgdir(RareEvents), "examples")
+include(joinpath(examples_dir, "standard_map.jl"))
+
+dt = 1.0
+alg_kwargs = ();
+FT = Float64
+K_sm =2.0
+model = StandardMap{Float64}(K_sm)
+evolve_wrapper(tspan) = (u0) -> evolve_mapping_system(model, u0, tspan)
+
+
+# Set up rare event algorithm
+τ = 10.0
+nensemble = 300
+u0 = [[rand()*2π, rand()*2π, 1e-5, 1e-5] for _ in 1:nensemble]
+kvals = [0.0,-0.05,-0.1]
+ϵ = 0.0001
+metric(y; K= K_sm) = (y[4]*y[3]*2.0*(1.0+K*cos(y[1])))/(y[3]^2.0+y[4]^2.0)
+t_end = 150.0
+t_start = 0.0
+tspan = (t_start, t_end)
+plot1 = plot()
+plot2 = plot()
+extract_I(u) = u[2]
+extract_θ(u) = u[1]
+colors = [:cyan, :purple, :red]
+ms = [1.0, 1.0, 1.0]
+labels = ["k=0", "k = -0.05", "k = -0.1"]
+# Run GKLT
+for j in 1:2
+    k = kvals[j]
+
+    sim = RareEventSampler{Float64}(dt, τ, tspan,u0, evolve_wrapper, nensemble, metric, k, ϵ, j)
+    run!(sim);
+    flatten_θ = extract_θ.(sim.ensemble[1])
+    flatten_I = extract_I.(sim.ensemble[1])
+    flatten_λ = -sum(metric.(sim.ensemble[1])*dt)/(tspan[2]-tspan[1])
+    for i in 2:nensemble
+        flatten_θ = vcat(flatten_θ, extract_θ.(sim.ensemble[i]))
+        flatten_I = vcat(flatten_I, extract_I.(sim.ensemble[i]))
+        flatten_λ = vcat(flatten_λ, sum(metric.(sim.ensemble[i])*dt)/(tspan[2]-tspan[1]))
+    end
+    scatter!(plot1, flatten_θ, flatten_I, color = colors[j], ms = ms[j], label = labels[j], markerstrokewidth=0)
+    plot!(plot2, flatten_λ, seriestype=:stephist, label = labels[j], bins = 10)
+end
+plot!(plot1, xlabel = "θ", ylabel = "I", xticks = ([0,π,2π],["0","π","2π"]), yticks = ([0, π, 2π],["0","π","2π"]))
+plot!(plot2, xlabel = "Finite time estimate of Lyapunov exponent λ", ylabel = "Relative Frequency")
+savefig(plot1, "standard_map.png")
+savefig(plot2, "lyapunov_histograms.png")