Computing velocity and rotation of rigid body after application of force

The motion of the body can be modeled by the system of first-order ordinary differential equations in vector/matrix: begin{align} J&frac{dvec{Omega}}{dt} , = , J,vec{Omega} times vec{Omega} , + , vec{X}_Ftimes vec{F}\ &\ &frac{d}{dt} U , = , U, big(vec{Omega} times cdotbig)\ &\ &frac{dvec{x}_G}{dt} , = ,vec{v}_G \ &\ &frac{dvec{v}_G}{dt} , = , U,vec{F} end{align} for the unknow vector and matrix functions $$vec{x}_G = vec{x}_G(t) in mathbb{R}^3,,,,,,vec{v}_G = vec{v}_G(t) in mathbb{R}^3, ,,,,, U = U(t) in text{SO}(3), ,,,,, vec{Omega} = vec{Omega}(t) in mathbb{R}^3$$ where the vector $vec{Omega}$ is the angular velocity of the body with respect to the body fixed frame (origin at the center of mass), $J$ is the 3 by 3 constant inertia matrix of the body, calculated in the body fixed frame and relative to the center of mass. The vector $vec{X}_F$ is the constant vector, in the body fixed frame, pointing from the center of mass $G$ of the body to the point of application of the force. The vector $vec{F}$ is the constant vector, also written in body fixed frame coordinates, of the force applied to a fixed point on the body. The change from body fixed coordinate frame to external inertial coordinate frame is given by the congruence transformation (rotation around center of mass plus translation along the vector $vec{x}_G$ in the inertial (external) coordinate system pointing from its origin to the center of mass of the body) $$vec{x} , = vec{x}_G , + , U, vec{X}$$ where the time dependent matrix function $U = U(t) in text{SO}(3)$ is an orthogonal matrix. The state of the body, also called the configuration of the body, is thus described by the pair $(,vec{x}_G,,, U,)$.
By $big(vec{Omega} times cdotbig)$ I mean the anti-symmetric 3 by 3 matrix with the property that for every three vector $vec{X}$ $$big(vec{Omega} times cdotbig), vec{X} , = , vec{Omega} times vec{X}$$

Now, the first vector equation is decoupled from the rest. There are many approaches how to integrate it, including geometric Lie-Poisson approximate integration or even in the case when the body has two equal axes of inertia it is exactly solvable with analytic formulas. But you can be just quick and dirty here and simply use explicit Runge-Kutta 4th order integrator with a step $frac{h}{2}$ giving you a sequence of angular velocities $vec{Omega}_t , , : ,, t = t_0 + nfrac{h}{2}$ with $n=0, 1, 2, ...$.

After that, when you have $vec{Omega}_t$ the second matrix equation is in fact decoupled from the third and forth, so maybe use as integration simple piece-wise rotations with step $frac{h}{2}$ to it as follows: if you have already calculated the rotation matrix $U_t$ then $$U_{t+frac{h}{2}} , = ,U_t , text{exp}left({frac{h}{2}big(vec{Omega}_t times cdotbig)}right) $$
where $text{exp}left({frac{h}{2}big(vec{Omega} times cdotbig)}right)$ is matrix exponentiation, which in the case of an anti-symmetric matrix $big(vec{Omega} times cdotbig)$ produces an orthogonal matrix. The exponentiation can be carried out explicitly by Rodrigues' formula. AS a result, you obtain a sequence of rotational (orthogonal) matrices $U_t , , : ,, t = t_0 + nfrac{h}{2}$ with $n=0, 1, 2, ...$.

Finally, apply the explicit Runge-Kutta 4th order to the last two vector equations because you know $U$ at time increments of $frac{h}{2}$.

This scheme will allow you to propagate the initial state of the body $big(vec{x}_G(0), ,, U(0)big)$ to the final one $big(vec{x}_G(t_1), ,, U(t_1)big)$.

The 3D motion of a rigid body requires some familiarity with vector and matrix operations.

1. The rotation about the center of mass does not have an analytical solution.

   Suppose at $t=0$ the body has orientation matrix $mathbf{R}_0$ as well as rotational velocity $boldsymbol{omega}_0$. This means the 3×3 mass moment of inertia tensor initially is calculated with the following matrix equation $mathbf{I}_0 = mathbf{R}_0 , mathbf{mathcal{I}}_{rm body} mathbf{R}_0^top$ where $mathbf{mathcal{I}}_{rm body}$ is the mass moment of inertia on body riding coordinates. Additionally, the angular momentum vector starts at $boldsymbol{H}_0 = mathbf{I}_0, boldsymbol{omega}_0$ if the initial rotational motion is $boldsymbol{omega}_0$.

   The body orientation matrix $mathbf{R}$ changes with time, according to the differential equation

   $$ tfrac{rm d}{{rm d}t} mathbf{R}(t) = boldsymbol{omega}(t) times mathbf{R}(t) tag{1}$$

   As the body continues to change orientation $mathbf{R}(t)$, if there was no external torque applied angular momentum vector $boldsymbol{H}$ would be conserved producing an expression for the rotation velocity of

   $$ begin{aligned}boldsymbol{omega}(t) & =mathbf{I}^{-1}boldsymbol{H}\ & =left(mathbf{R}(t),mathbf{mathcal{I}}_{{rm body}},mathbf{R}^{top}(t)right)^{-1}boldsymbol{H}_{0}\ & =left(mathbf{R}(t),mathbf{mathcal{I}}_{{rm body}}^{-1},mathbf{R}^{top}(t)right),boldsymbol{H}_{0}\ & =left(mathbf{R}(t),mathbf{mathcal{I}}_{{rm body}}^{-1},mathbf{R}^{top}(t)right),underbrace{mathbf{R}_{0},mathbf{mathcal{I}}_{{rm body}},mathbf{R}_{0}^{top}boldsymbol{omega}_{0}}_text{constant} end{aligned} tag{2}$$

   The expression (1) does not have an analytical solution for $boldsymbol{R}(t)$ in light of the torque-free equation (2).

   But since the torque applied about the center of mass is variable, $boldsymbol{tau}(t) = mathbf{R}(t) left( boldsymbol{mathcal{c}} times boldsymbol{mathcal{F}} right)$ where $boldsymbol{mathcal{c}}$ is the location of the force in body coordinates, and $boldsymbol{mathcal{F}}$ the force in body coordinates. Now this torque modifies the angular momentum vector by the following differential equation

   $$ tfrac{rm d}{{rm d}t} boldsymbol{H}(t) = boldsymbol{tau}(t) tag{3} $$

   As a result (2) becomes a lot more complex with

   $$ boldsymbol{omega}(t) = left( mathbf{R}(t) mathbf{mathcal{I}}_{rm body}^{-1} mathbf{R}^top(t) right) boldsymbol{H}(t) tag{4}$$

2. The motion of the center of mass depends on the orientation of the body.

   The force aligned with the body is $mathcal{F}$, and in world coordinates $boldsymbol{F}(t) = mathbf{R}(t),boldsymbol{mathcal{F}}$ and from Newton's second law the acceleration of the center of mass is

   $$boldsymbol{a}(t) = tfrac{1}{m}, boldsymbol{F}(t) = tfrac{1}{m},mathbf{R}(t) boldsymbol{mathcal{F}} tag{5} $$

   And the position and velocity calculated from the following differential equations

   $$begin{aligned} tfrac{rm d}{{rm d}t} boldsymbol{v}(t) &= boldsymbol{a}(t) \ tfrac{rm d}{{rm d}t} boldsymbol{r}(t) &= boldsymbol{v}(t) end{aligned} tag{6}$$

3. Special Cases

   After extensive research NASA in the 60's developed some analytical solutions of the above equations for the special case of a nearly axisymmetric body, which is quite a feat considering the complexity of the equations above.

Well i would assume you need the velocity of a particular particle in the system.

So the best way to proceed is-

First find the velocity of the centre of mass. This can be found directly by applying Newton's 2nd law for the applied force.

Then you can find the velocity of the particle in the centre of mass frame. This is done by finding torque due to the force, and using that in the adopted form of Newton's 2nd law to find angular velocity. With this you can find the velocity of the particle w.r.t centre of mass.

And finally you can use these 2 relations to find velocity of the point in the ground frame