function dynamics!(
	state,
	prev_state,
	params,
	states,
	plot_data,
)
	# Unpack the parameters
	N, K, beta, A, dt, num_steps = params

	for j in 1:num_steps
		next_state = zeros(N)

		# update the left particle (set to rest)
		state[1] = 0

		# update the right particle (set to rest)
		state[N] = 0

		# update the middle particles
		for i in 2:N-1
			a = 2 * state[i] - prev_state[i]
			b = dt * dt * K * (state[i+1] - 2 * state[i] + state[i-1])
			c = dt * dt * beta * ((state[i+1] - state[i])^3 + (state[i-1] - state[i])^3)

			this_mass = 1
			if i == Int(N / 2)
				# time = i * dt
				# push!(plot_data, (state[i]))
				this_mass = 1
			end
			next_state[i] = (a + b + c) / this_mass
		end

		# push!(states, next_state)

		# update the previous state
		for i in 1:N
			prev_state[i] = state[i]
		end
		# update the current state
		for i in 1:N
			state[i] = next_state[i]
		end

		# if the middle mass is close to 2, print the time and position
		if state[Int(N / 2)] > 1.995
			println("Time: ", j * dt, " Positions: ", state[Int(N / 2)])
			println("Other Positions: ", state, "\n")
			push!(states, (j * dt, state))
		end

		if (j * dt % 100000) == 0
			println("TIME UPDATE: ", j * dt)
		end
	end
end

function get_initial_state(
	N,
	modes,
	beta,
	A,
)
	state = zeros(N)
	# amp = A * sqrt(2 / (N - 1))
	# for i in 2:N-1
	# 	state[i] = amp * sin((modes * π * (i - 1)) / (N - 1))
	# end
	# set the middle to be the largest
	state[Int(N / 2)] = 2
	return state
end

function run_simulation(
	N,
	modes,
	beta,
	A,
	dt,
	num_steps,
	plot_data,
)
	states = []
	state = get_initial_state(N, 1, beta, A)
	prev_state = zeros(N)
	for i in 1:N
		prev_state[i] = state[i]
	end
	params = (N, modes, beta, A, dt, num_steps)
	dynamics!(state, prev_state, params, states, plot_data)
	return states
end

# Run the simulation
N = 10 # number of masses
beta = 0 # cubic string spring
K = 100 # spring constant
A = 10 # amplitude
modes = 3 # number of modes to plot
final_time = 10000000 # seconds
dt = 0.001 # seconds
num_steps = Int(final_time / dt)
params = (N, K, beta, A, dt, num_steps)
plot_data = []
s = run_simulation(N, K, beta, A, dt, num_steps, plot_data)
println("\nStates of interest:", s)

# plot the inital positions and the final positions
# using Plots
# plot(get_initial_state(N, 1, beta, A), label = "Initial", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", title = "First Three Modes")
# plot!(s[end], label = "Final", marker = :circle)

# plot the middle of the lattice over time
# time_steps = [i * dt for i in 1:num_steps]
# plot(time_steps, plot_data, label = "Middle Mass Over Time", xlabel = "Time", ylabel = "Displacement")
# savefig("t/plot_data.png")

# print the points close to the origional displacement
# output = []
# for i in 1:length(plot_data)
# 	if plot_data[i] > 1.975
# 		push!(output, i)
# 		println("Time: ", i * dt, " Position: ", plot_data[i])
# 	end
# end


# animate the s array of positions
# anim = @animate for i in 1:length(s)
# 	plot(s[i], label = "t = $(round(i * dt, digits = 3))", marker = :circle, xlabel = "Mass Number", ylabel = "Displacement", ylim = (-3, 3))
# end
# mp4(anim, "t/plz.mp4", fps = 30)
# println("Output: ", output)