Single harmonic oscillator has been
ubiquitous in both classical as well as quantum mechanics. The harmonic
oscillator is governed by the potential:
$U(x) = \dfrac{1}{2}kx^2 $
It is quite easy to see that the equations of motion for this case becomes:
$\dot(x) =p$
$\dot(p) =-x$
Any deviation in potential of the first equation causes the oscillator to forgo its "harmonicity" and thereby, enters anharmonic nature. Anharmonicity plays an important role in several toy models of thermal conduction like Fermi-Pasta-Ulam and $\phi^4$ chain. The focus of this short write-up is to show how the numerical solution changes when the system is no longer harmonic. For simplicity we assume that the anharmonic potential is of the form:
$U(x) = \dfrac{1}{2}x^2 + \dfrac{1}{3}x^3 + \dfrac{1}{4}x^4 +...+ \dfrac{1}{2n+1}x^{2n+1} $
Note, we terminate at the index 2n+1 because otherwise the solutions
diverge. The logic of such divergence is easy to understand: The force
is always repulsive and further the oscillator moves, greater is the
repulsive force. The MATLAB codes are attached below for the interested
readers. The four resulting solutions are given below:
The MATLAB codes are given below:



