Accurate 2D collision simulation in Rust

I have created a 2D physics simulation in rust that works for the most part. It involves creating balls which collide with each other and their container. I have looked at a lot of resources and found that this and this have been the most helpful (there are suprisingly few resources for collisions of two moving balls in 2D).

The maths I used to model the collision is specified in this pdf: https://www.vobarian.com/collisions/2dcollisions2.pdf.

Below is a video of what problem occurs:

As you can see, some of the balls will get stuck inside of each other, rotate a little and then pop out of each other. It is very unsual behaviour and shouldn't be possible given that I have (as far as I am aware) implemented the maths correctly. Below is my collision function, which is iterated through on every ball.

fn collide(ball_1: &mut Ball, ball_2: &mut Ball) {
// Checking for collisions in next frame to stop balls from overlapping, so that they won't get "trapped"
let next_pos_1 = ball_1.pos + ball_1.vel;
let next_pos_2 = ball_2.pos + ball_2.vel;

let distance =((next_pos_1[0]-next_pos_2[0]).powf(2.0)+(next_pos_1[1]-next_pos_2[1]).powf(2.0)).sqrt(); 
if distance < ball_1.radius + ball_2.radius {
    let normal_vector = ball_2.pos - ball_1.pos;
    
    let unit_normal = normal_vector / (normal_vector[0].powi(2) + normal_vector[1].powi(2)).powf(0.5);
    let unit_tangent = Vec2::new(-unit_normal[1], unit_normal[0]);
    
    let v_1_in_norm_dir = unit_normal.dot(ball_1.vel);
    let v_2_in_norm_dir = unit_normal.dot(ball_2.vel);
    let v_1_in_tan_dir = unit_tangent.dot(ball_1.vel);
    let v_2_in_tan_dir = unit_tangent.dot(ball_2.vel);

    let new_v_1_in_norm_dir = v_2_in_norm_dir;
    let new_v_2_in_norm_dir = v_1_in_norm_dir;

    ball_1.vel = new_v_1_in_norm_dir * unit_normal + v_1_in_tan_dir * unit_tangent;
    ball_2.vel = new_v_2_in_norm_dir * unit_normal + v_2_in_tan_dir * unit_tangent;

    let next_pos_1 = ball_1.pos + ball_1.vel;
    let next_pos_2 = ball_2.pos + ball_2.vel;
    let distance =((next_pos_1[0]-next_pos_2[0]).powf(2.0)+(next_pos_1[1]-next_pos_2[1]).powf(2.0)).sqrt(); 

    // Edge case where the ball's collision will cause one to push the other, so just push them back apart
    if distance >= ball_1.radius +ball_2.radius {
        let overlap = ((ball_1.radius + ball_2.radius) - distance)/2.0;

        let m = (ball_1.pos[1]-ball_2.pos[1])/(ball_1.pos[0]-ball_2.pos[0]);
        let theta = m.atan().abs();
        let y_change = overlap*theta.sin();
        let x_change = overlap*theta.cos();
        
        if ball_1.pos[0] <= ball_2.pos[0] {
            ball_1.pos[0] -= x_change;
            ball_2.pos[0] += x_change;
            ball_1.pos[1] -= y_change;
            ball_2.pos[1] += y_change;
        } else {
            ball_1.pos[0] += x_change;
            ball_2.pos[0] -= x_change;
            ball_1.pos[1] += y_change;
            ball_2.pos[1] -= y_change;
        }
    }
    
}   

}

The last if statement is trying to combat my current problem. My current theory on what is occuring is that sometimes, when the velocities of the balls are calculated correctly, these velocities will actually force the balls into one another (because the velocity of one ball points towards the other ball). When the velocities are then added onto the positions at the end of the function, it causes one ball to overlap with the other. This is problematic as this function is supposed to only trigger when the balls are about to overlap or are already overlapping. This is therefore causing the balls to overlap more and more and just create a kind of loop, which eventually ends. It is causing the collision function to do the opposite of what it is supposed to do.

I might be wrong about this theory but it is quite difficult to debug. I am not a physicist (clearly) and so I am struggling with figuring out what the best thing to do here is. I want to keep the accurate collisions and processing speed but lose this glitchy behaviour.

I couldn't find any posts talking about this but please feel free to give me any advice as I am feeling a bit lost! Thanks :)

I have mostly figured it out after spending a long time debugging. It appears that there were some problems with the solution for resolving all collisions each frame. I was using a grid solution where each ball only checks if it has collided with balls in its grid. But this was creating problems and so now a simpler solution has been implemented. The whole code can be seen below:

use macroquad::prelude::*;
use ::glam::Vec2;

const G: f32 = 0.0; // Gravity constant
const BG_COLOR: Color = BLACK;
const GRID_SIZE: [i32; 2] = [10, 10];

struct Ball {
    color: Color,
    radius: f32,
    pos: Vec2,
    vel: Vec2,
    acc: Vec2,
    boundary_cor: f32, // COR that the ball has with boundaries  https://en.wikipedia.org/wiki/Coefficient_of_restitution
}

impl Ball {
    fn update(&mut self, screen: [f32; 2]) {
        self.vel += self.acc;
        self.pos += self.vel;

        // collision w boundary detection
        if self.pos[0]+self.radius > screen[0] {
            self.vel[0] = -self.vel[0]*self.boundary_cor;
            self.pos[0] = screen[0]-self.radius;
        } else if self.pos[0]-self.radius < 0.0 {
            self.vel[0] = -self.vel[0]*self.boundary_cor;
            self.pos[0] = self.radius+0.1;
        }
        if self.pos[1]+self.radius > screen[1]{
            self.vel[1] = -self.vel[1]*self.boundary_cor;
            self.pos[1] = screen[1]-self.radius;
        } else if self.pos[1]-self.radius < 0.0 {
            self.vel[1] = -self.vel[1]*self.boundary_cor;
            self.pos[1] = self.radius+0.1;
        }
        //collision w cursor detection
        let cursor = mouse_position();
        if cursor.0 > self.pos[0]-self.radius && cursor.0 < self.pos[0]+self.radius {
            if cursor.1 > self.pos[1]-self.radius && cursor.1 < self.pos[1]+self.radius {
                self.vel[1] = -self.vel[1];
                self.vel[0] = -self.vel[0];
            }
        }
    }
}

fn collide(ball_1: &mut Ball, ball_2: &mut Ball) {
    // It's okay to check current frame because balls have just been updated but not displayed, so can still stop them from being inside each other for this frame

    let distance = ((ball_1.pos[0]-ball_2.pos[0]).powf(2.0)+(ball_1.pos[1]-ball_2.pos[1]).powf(2.0)).sqrt(); 
    if distance <= ball_1.radius + ball_2.radius {
        

        // Resolve vels
        let v_1 = ball_1.vel;
        let v_2 = ball_2.vel;

        ball_1.vel += (v_2 - v_1).dot(ball_2.pos - ball_1.pos) / distance / distance * (ball_2.pos - ball_1.pos);
        ball_2.vel += (v_1 - v_2).dot(ball_1.pos - ball_2.pos) / distance / distance * (ball_1.pos - ball_2.pos);

        // Removes overlap
        // TODO: Make sure this overlap removal doesn't push one of the balls a bit off screen, can cause problems
        let overlap = (ball_1.radius + ball_2.radius) - distance;
        let dir = (ball_1.pos - ball_2.pos).clamp_length(overlap/2.0, overlap/2.0);
        ball_1.pos += dir;
        ball_2.pos -= dir;
    }  
}


#[macroquad::main("WAZZA")]
async fn main() {
    let color_choices = [RED, BLUE, YELLOW, ORANGE];
    let mut balls: Vec<Ball> = Vec::new();
    for i in 0..5 {
        for j in 0..5 {
            balls.push(Ball {
                color: color_choices[i%4],
                radius: 20.0,
                pos: Vec2::new((i as f32)*((screen_width() as f32)/10.0), (j as f32)*((screen_height() as f32)/10.0)),
                vel: Vec2::new(i as f32, j as f32),
                acc: Vec2::new(0.0, G),
                boundary_cor: 0.9,
            });
        }
    }
    loop {
        clear_background(BG_COLOR); //Screen is cleared whether or not function is called so no performance reduction
        let screen = [screen_width(), screen_height()];
        
        for ball in 0..balls.len() { // Physics update all balls
            balls[ball].update(screen);
        }

        // Ball to ball collision detection
        for ball in 0..balls.len() {
            for other_ball in ball+1..balls.len() {
                let (left, right) = balls.split_at_mut(other_ball);
                collide(&mut left[ball], &mut right[0]);
            }
        }

        let mut total_momentum = 0.0;
        for ball in 0..balls.len() {
            draw_circle(balls[ball].pos[0], balls[ball].pos[1], balls[ball].radius, balls[ball].color);

            total_momentum += balls[ball].vel.length();
        }
        println!("{}", total_momentum);

        draw_fps();
        next_frame().await
    }
}