I've been trying to wrap my head around this problem for some time and I just can't seem to get it straight.

What I am trying to achieve:
I want to program a recursive model that calculates the position of a moving body in a fluid after a certain time.

What i have done already
I got as far as to simulate the movement due to a force applied at the center of the object.
In pseudo-code

Spoiler

for the speed at t = n
//calculate the drag force due to v_{n-1}Wikipedia reference
F_{d} = 0.5*ρ*v_{n-1}^{2}*C_{d}*A //
// add any external force that may have been applied when n-1<=t<n and the gravitational force
F_{tot}= F_{d}+F_{ext}+F_{g}
// calculate the acceleration of the object at n
a_{n}=F_{tot}/m
// assume the acceleration constant for dt so the velocity becomes
v = v_{n-1} + a_{n}*dt
// calculate the new coordinates of the object based on v_{n}

This works fine as long as the time step dt is small enough (0<dt<0.00001)
The amount of calculations is too high. So i decided to approximate some more by limiting the rate of change of the drag force F_{d.est} = MAX(const,F_{d})
This leads to a decent approximation with about a factor 100 less calculations per second (0<dt<0.001), which i think is sufficient. I observed a 10% disparity between the original and the approximation after 10^7 iterations.

What i'm struggling with
1) I'm stumped on how to model the angular motion when a force is applied a distance d from the center of the object
2) Is there a way to make the approximate the whole model more precisely but with the same dt = 0.001

If anyone has tips on how I can improve I would greatly appreciate it.

For off-center forces, you will have to compute the produced torques and solve another equation for angular motion.

To compute things more efficiently, you will need to learn some numerical methods for ordinary differential equations. Your method, by the way, is known as Euler's method.

My dynamics are not what they ought to be. Could you provide an example calculation for the torque and the angular motion equations?
Is there a numerical method for ODE s that I can readily apply or if not, is there free literature I can read up in?

First, you need to determine the center of mass and the moment(s) of inertia of your body. Then you need to determine the net torque of the forces (relative the center of mass). Then you will have an equation very similar to Newton's second law. Look up the details on Wikipedia, if unclear, ask for help here.

To solve an ODE, I would try by default one of the Runge-Kutta integrators. Wikipedia is your friend again.

So for the off-center forces:
// calculate all moments due to forces that apply on the body (drag, external, gravity = 0 due to symmetry)
F_{d}= 0.5*A*C_{d}*ρ*ω_{n-1}^{2}
M_{drag} = F_{d}*r where r = radius of the object(assuming the object is symmetrical)
M_{ext} = F_{ext}*x where x is the distance to to centre of mass
// calculate angular acceleration
α = M_{tot}/I
// calculate angular speed
ω_{n}=ω_{n-1}+α*dt
// calculate new angle of the body
θ_{n}=θ_{n-1}+ω_{n}*dt

For the runge kutta methods. If I understood correctly I would want to calculate the (angular) velocity with it by doing above calculations 4 times with
k_{1} : euler
k_{2} : dt = dt/2 and v_{n-1} = v_{n-1} + k1*dt/2
k_{3} : dt = dt/2 and v_{n-1} = v_{n-1} + k2*dt/2
k_{4} : dt = dt and v_{n-1} = v_{n-1} + k3*dt

and then v_{n} = v_{n-1} + dt/6(k_{1}+2k_{2}+2k_{3}+k_{4})

and the new coordinates are then p_{n} = p_{n-1} + v_{n}*dt

The first part of that looks OK, but the last part
p_{n} = p_{n-1} + v_{n}*dt

is throwing all improved accuracy away, because it is just the same as Euler's forward difference method!

You need to rewrite your second-order differential equation as two coupled first order equations and then solve both of them together using R-K. Any examples of solving equations of motion usiing R-K (on the web or in a textbook) should show you how to do it right. (Or search PF for other threads about RK, this sort of question comes up fairly often).

Also, it's not obvious from your equations that you will always have the correct sign for the drag force. Remember the force changes direction if the velocity is reversed, but ##v^2## does mot change sign when ##v## goes from positive to negative. Maybe you want to calcule ##v | v | ## instead (i.e. ##v## times the absolute value of ##v##.)