well seems to be related to this thread. Note that Levenberg Marquardt is not the robustest Method for such kind of problems (multiple local minima).
Optim.jl provides a simple SA-Optimizer which tend to be more robust for this task
using Optim
loss(u) = sum((data .- model(t,u)).^2)
res = optimize(loss, p0_bounds, SimulatedAnnealing())
# Plot result
scatter(t, data)
plot!(t, model(t, Optim.minimizer(res)))
Or you can try the awesome BlackBoxoptim.jl
EDIT: To elaborate a bit on this:
By simply plotting the value of loss(u) around the p₀ you can visualize how challenging your problem actually is:
you can see multiple local minima next to the global one