elastic collisions of two disks and the surrounding walls in matlab

main.m
% Parameters
m1 = 1; % mass of disk 1
m2 = 1; % mass of disk 2
v1_initial = [1, 0]; % initial velocity of disk 1 [vx, vy]
v2_initial = [-1, 0]; % initial velocity of disk 2 [vx, vy]
r1 = 0.5; % radius of disk 1
r2 = 0.5; % radius of disk 2
e1 = 1; % coefficient of restitution for disk 1 (perfectly elastic)
e2 = 1; % coefficient of restitution for disk 2 (perfectly elastic)
width = 10; % width of the walls
height = 10; % height of the walls
timestep = 0.01; % timestep for the simulation
simulation_duration = 10; % duration of the simulation

% Initial positions
p1 = [0, 0]; % position of disk 1 [x, y]
p2 = [5, 0]; % position of disk 2 [x, y]

% Simulation loop
t = 0; % current time
while t < simulation_duration
    % Calculate new velocities after collision
    v1_initial_x = v1_initial(1);
    v1_initial_y = v1_initial(2);
    v2_initial_x = v2_initial(1);
    v2_initial_y = v2_initial(2);
    v1_final_x = ((m1 - m2) * v1_initial_x + 2 * m2 * v2_initial_x) / (m1 + m2);
    v1_final_y = ((m1 - m2) * v1_initial_y + 2 * m2 * v2_initial_y) / (m1 + m2);
    v2_final_x = ((m2 - m1) * v2_initial_x + 2 * m1 * v1_initial_x) / (m1 + m2);
    v2_final_y = ((m2 - m1) * v2_initial_y + 2 * m1 * v1_initial_y) / (m1 + m2);
    
    % Update positions of the disks
    p1 = p1 + v1_initial * timestep;
    p2 = p2 + v2_initial * timestep;
    
    % Check for collision with walls
    if p1(1) < -width/2 + r1 || p1(1) > width/2 - r1
        v1_final_x = -e1 * v1_final_x;
    end
    if p1(2) < -height/2 + r1 || p1(2) > height/2 - r1
        v1_final_y = -e1 * v1_final_y;
    end
    
    if p2(1) < -width/2 + r2 || p2(1) > width/2 - r2
        v2_final_x = -e2 * v2_final_x;
    end
    if p2(2) < -height/2 + r2 || p2(2) > height/2 - r2
        v2_final_y = -e2 * v2_final_y;
    end

    % Update velocities for next iteration
    v1_initial = [v1_final_x, v1_final_y];
    v2_initial = [v2_final_x, v2_final_y];
    
    % Increment time
    t = t + timestep;

    % Plot current position of the disks
    hold off;
    plot(p1(1), p1(2), 'ro', 'MarkerSize', r1*100);
    hold on;
    plot(p2(1), p2(2), 'bo', 'MarkerSize', r2*100);
    rectangle('Position', [-width/2, -height/2, width, height]);
    axis([-width/2, width/2, -height/2, height]);
    axis equal;
    drawnow;
end

2293 chars
69 lines