>The model I am using for QBO is that there is a characteristic resonant frequency associated with the wave equation, but that an exact biannual modulation is applied.
In the code it seems you are not using the above equation but the following Hill equation
$$y''(x) + (const1+ const2*cos(4 \pi x + const3) + const4*cos(0.5195 \pi x + const5))*y(x)=0$$
I am not sure (anymore) how the solutions of Hills equation look like, but I would think that if the solution is periodic then it would have periodicities as in the given fourier expansion, so this points to a periodicity
of 1/0.5195 and not 1/0.5 (which would be the exact biannuity).
I assume that "First" plots out y(x) of NDSolve in particular I don't have Mathematica to check. Moreover in your plot you have that strange $y(x+const*sin(const x +....)
is that what you call a filter?
Where do you have the QBO data from? Did you compare that with the one in the blogpost you had linked to?