function ukk = rk4_feval(f,tk,uk,omega,gamma,Dt)

k1 = feval(f,tk,uk,omega,gamma);
k2 = feval(f,tk,uk+0.5*Dt*k1,omega,gamma);
k3 = feval(f,tk,uk+0.5*Dt*k2,omega,gamma);
k4 = feval(f,tk,uk+Dt*k3,omega,gamma);

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