Here's a quick way of doing it with one matrix multiplication, you just need to precompute cos(ω) and sin(ω) with a bit of intuition:

   [ x(n + 1) ] = [cos(ω), -sin(ω)] [ x(n) ]
   [ y(n + 1) ]   [sin(ω),  cos(ω)] [ y(n) ]

   where 

   x(n) = cos(ω n)
   y(n) = sin(ω n)
   x(0) = 1
   y(0) = 0
   ω = frequency (radians) * time_step
   

There are two ways to look at this, from the systems perspective this is computing the impulse response of a critically stable filter with poles at e^(+/-jω). Geometrically, it's equivalent to starting at (1, 0) and rotating around the unit circle by ω radians each time step. You can offset to arbitrary phases.

It's not suitable for numerous cases (notably if the frequency needs to change in real time) but provided that the sum of the coefficients sum to `1.0` it's guaranteed stable and will not alias.