#### Runge Kutta Integration

P
"""
runge_kutta_integration(f::Function, x0::Real, y0::Real, h::Real, x_stop::Real)

Numerically solve initial value problems of the form y' = f(x, y) to find y(x) using the Runge-Kutta 4th order integration scheme.

Starting with the differential equation \\frac{\\mathrm{d}y}{\\mathrm{d}x} = y' = f(x, y) and the initial condition y(x_0) = y_0, each step calculates 4 approximations of the gradient
math
\\begin{align*}
k_1 &= f(x_n, y_n),\\\\
k_2 &= f(x_n + \\frac{h}{2}, y_n + k_1\\frac{h}{2}),\\\\
k_3 &= f(x_n + \\frac{h}{2}, y_n + k_2\\frac{h}{2}),\\\\
k_4 &= f(x_n + h, y_n + k_3h),
\\end{align*}

and uses the weighted average of them,

math
\\bar{k} = \\frac{k_1 + 2k_2 + 2k_3 + k_4}{6},


to find the approximate value of y(x_{n+1}) and update x and y accordingly

math
\\begin{align*}
x &\\rightarrow x_{n+1} = x_n + h\\\\
y &\\rightarrow y_{n+1} = y_n + h\\bar{k}.
\\end{align*}


# Arguments:
- f: The function y' = f(x, y) to solve for y(x).
- x0: The starting value of x.
- y0: The starting value of y.
- h: The step size, should be >0.
- x_stop: The final value of x to solve to, i.e. integrate over the range [x0, x_stop]

# Examples
julia
julia> # If you have a constant slope of y' = 1, the analytic solution is y = x
julia> runge_kutta_integration((x, y)->1, 0, 0, 1, 3)
([0.0, 1.0, 2.0, 3.0], [0.0, 1.0, 2.0, 3.0])

julia> # Consider solving y' = exp(x), which has the analytic solution y = exp(x).
julia> runge_kutta_integration((x, y)->exp(x), 0, 1., 0.1, 0.1)
([0.0, 0.1], [1.0, 1.105170921726329])

julia> exp.([0.0, 0.1])
2-element Vector{Float64}:
1.0
1.1051709180756477



# References
- [https://en.wikipedia.org/wiki/Initial_value_problem](https://en.wikipedia.org/wiki/Initial_value_problem)
- [https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods](https://en.wikipedia.org/wiki/Initial_value_problem)

# Contributors:
- [E-W-Jones](https://github.com/E-W-Jones)
"""
function runge_kutta_integration(
f::Function,
x0::Real,
y0::Real,
h::Real,
x_stop::Real,
)
h > 0 || throw(DomainError(h, "The step size h should be >0."))

x = Float64(x0)
y = Float64(y0)
output_x = [x]
output_y = [y]

while x < x_stop
k1 = f(x, y)
k2 = f(x + h / 2, y + k1 * h / 2)
k3 = f(x + h / 2, y + k2 * h / 2)
k4 = f(x + h, y + k3 * h)

y += h * (k1 + 2k2 + 2k3 + k4) / 6
x += h

push!(output_x, x)
push!(output_y, y)
end

return output_x, output_y
end