function ukk = m_theta(theta,omega,gamma,Dt,uk)

ukk = 1/(1+2*gamma*theta*Dt+(theta*Dt*omega)^2)*...
	 [1+2*gamma*theta*Dt, theta*Dt; -theta*Dt*omega^2, 1]*...
	 [1, (1-theta)*Dt; -(1-theta)*Dt*omega^2, 1-2*gamma*(1-theta)*Dt]*uk;