Numerical method for harmonic oscillator with jumping constant

Let $$k_1 \neq k_2$$ be positive reals, $$t_0 > 0$$ and consider the following Cauchy problem in $$[0,+\infty)$$: $$\begin{cases} y(t) + k(t)y''(t) = 0 \newline y(0) = 1/\sqrt{k_1} \newline y'(0) = 0, \end{cases}$$ where $$$$k(t) = \begin{cases} k_1 \hspace{1 cm} 0 \leq t \leq t_0 \newline k_2 \hspace{1 cm} t > t_0 \end{cases}$$$$

It is clear that there exists a unique $$C^1$$ solution $$y$$ in $$[0,+\infty)$$, obtained by gluing two smooth solutions of the equation in $$[0,t_0]$$ and $$[t_0,+\infty)$$. An explicit computation gives: $$$$y(t) = \begin{cases} \dfrac{1}{\sqrt{k_1}} \mathrm{cos} \bigg(\dfrac{t}{\sqrt{k_1}} \bigg) \hspace{1 cm} 0 \leq t \leq t_0; \newline A \ \mathrm{cos}\bigg( \dfrac{t}{\sqrt{k_2}} + \phi \bigg) \hspace{1 cm} t > t_0, \end{cases}$$$$ where $$A$$ and $$\phi$$ are obtained imposing the $$C^1$$ condition. With further computations one discovers that the amplitude is constant, i.e. $$A = 1/\sqrt{k_1}$$, if and only if $$t_0 = n \pi \sqrt{k_1}$$ for some $$n \in \mathbb N$$.

I am searching for a good numerical method for the problem in the constant-amplitude case. I have tried several classical schemes, but I discovered that when $$k_2 >> k_1$$ the second amplitude also grows bigger. Any suggestion?

• Try a symplectic integrator from DifferentialEquations.jl in Julia: docs.sciml.ai/dev then use events to handle the jump. You can also specify that you want to integration steps to include $t_0$ so you hit the jump dead-on Apr 24 '20 at 3:29
• @whpowell96 thank you, I'll have a look. Apr 24 '20 at 6:22
• @ChrisRackauckas Apr 24 '20 at 18:06