How to Calculate and Plot Fit Residuals from Numerical Integration to High Precision

The @. command is a macro that converts all the subsequent commands to broadcasted form:

help?> @.
  @. expr

  Convert every function call or operator in expr into a "dot call" (e.g.
  convert f(x) to f.(x)), and convert every assignment in expr to a "dot
  assignment" (e.g. convert += to .+=).

You could equivalently write that line as err[err .== 0] .= NaN, but I prefer to use the dot macro for conditional replacement, since the conditions can get pretty convoluted. You can read more about the @. macro in this blog post: More Dots: Syntactic Loop Fusion in Julia

The ‘problem’ is that your example is too easy for DiffEq.jl to integrate, so it needs very few timesteps to find an acceptable solution. Artificially clamping the timestep won’t really give you a clearer picture of the solver’s accuracy, since you wouldn’t choose such a fine discretization in day-to-day use. I would instead use a (smooth) line plot for the interpolated solution’s residual, overlaid with a scatter plot of the raw residual. Continuing on from my previous code,

t_fine = LinRange(tspan..., 1000)
u_DiffEq_fine = [solSpring1D(t)[2] for t in t_fine]
u_analytical_fine = @. 5cos(3t_fine)
err_fine = abs.(u_analytical_fine - u_DiffEq_fine)
@. err_fine[err_fine == 0] = NaN

plot(t_fine, err_fine, label="Interpolated solution")
scatter!(solSpring1D.t, err, label="Raw solution", 
     xlabel = "time", ylabel="Absolute error", yaxis=:log10, legend=:bottomright)

image