Nature of Code: Drag Forces

This might be more a physics question then a programming question but because it is described in “The Nature of Code” somebody might have an answer as well.
I am currently working on something, where I want to implement drag to a moving/falling object.
Following, what is written in TNoC (and what I think makes sense), I used the drag equation calculation a force and then applied that to my moving object.
I apply the constant force from gravity as well.

What I would expect is that the velocity is increasing to an equilibrium/constant value, where drag and gravity equal out
What happens is that the velocity is fluctuating around 0, quickly escalating to endless values.

Just a simple example, why:
Speed = 0 (Gravity +10, Drag -0)
Speed = 10 (Gravity +10, Drag -100)
Speed = -80 (Gravity +10, Drag +6400)
Speed = 6330 etc.

Because drag goes with the speed square it quickly dominates the forces defining acceleration. Obviously that should not happen, but I am not sure, what is wrong, because it is simply applying the physical equations.

I probably could limit the drag or use a drag constant so small, that it has a smaller effect, but I would like to use the exact equation, not just manipulating it to make it work.
I guess I miss something in the application of the drag equation?

 void force(){
   PVector gforce = new PVector(0,3);
   applyForce(gforce);
   PVector drag = velocity.copy();
   float cd = 1.2;
   float speed = velocity.mag();
   float dragmagnitude = cd * sq(speed);
   //drag.setMag(dragmagnitude);
   drag.mult(-1);
   drag.normalize();
   drag.mult(dragmagnitude);
   applyForce(drag);
  }
   
   void applyForce(PVector force_){
     PVector f = PVector.div(force_,mass);
     acceleration.add(f);     
   }
   
   void update(){
     velocity.add(acceleration);
     location.add(velocity);
     acceleration.mult(0);
   }

Yeah, you need a constant drag coefficient in there.

Normally,
force_drag = coeff_drag * fluid_density * area * 1/2 * speed^2

But this part: coeff_drag * fluid_density * area * 1/2

can all be bundled up into one constant. That you are getting a drag value that’s TEN TIMES more than the acceleration from gravity seems… excessive. Maybe try a cd value of of 0.01 instead.

I used cd as drag coefficient, which is usually 1.2 for plates.
Area and water density are also larger then 1, especially for larger objects, so using an extremely small constant, as you suggest, would only make sense for tiny objects.

I think the issue is, that the acceleration is just added to the velocity. That only works, if acceleration is significantly smaller than the velocity in absolute value.

Acceleration can be larger then velocity in value but then it acts faster and therefore has to be applied in smaller timesteps/portions.
So I would have to add it in tenth or hundredth of the total acceleration per frame.

Alternative would probably be calculating the equilibrium value.