function ukk = rk4(omega,gamma,Dt,uk)

A = [0, 1; -omega^2, -2*gamma];

k1 = A*uk;
k2 = A*(uk+0.5*Dt*k1);
k3 = A*(uk+0.5*Dt*k2);
k4 = A*(uk+Dt*k3);

ukk = uk + Dt*(k1/6+k2/3+k3/3+k4/6);