Projectile motion with drag. My problem is in the loop as the number of elements are not equal apparently on the left and right. Id appreciate detailed input ive been lost on this for a couple of days now

function Yoted 
clear
clc
format long g;
workspace
format compact 
fontSize = 20;

g = 9.81;                           % acceleration due to gravity in Y
y0 = 0;                             % intial position y
vHit = 50;                          % initial velocity magnitude
v0Max = vHit/10 + vHit;
v0Min = vHit - vHit/10;
v0 = randi([v0Min v0Max],1)         %randomly generated v0        
wHit = 100;                         % Angular Velocity (rad/s)
wMin = wHit - wHit/10; 
wMax = wHit + wHit/10;
Angular_w = randi([wMin wMax],1)    % randomly generate w
Angle = 30;                         % angle of strike
aMin = Angle - Angle/10;
aMax = Angle + Angle/10;
Alpha = randi([aMin aMax],1)        % Angle Alpha hit
c = 4.5*10^(-5);                    % Coefficient of drag
r = 21.5*10^(-3);                   % Radius of gold ball
m = 46*10^(-3);                     % Mass of golf ball 
rho_air = 1.2;                      % Density of air 


p = [0 0 0];
dt=0.1; 
t=[0.1:dt:200];
v = v0;                                    % Velocity Magnitude 

Velx = v(1)*cosd(Alpha)
Vely = v(1)*sind(Alpha)- g*t(1)
Vel(1,:) = [Velx Vely 0]                          % Velocity Vector 


FdragX(1) = v(1).*Vel(1).*c;                  % Drag force in X direction
FdragY(1) = v(1).*Vel(2).*c;                 % Drag force in Y direction
AccDragX(1) = -FdragX(1)./ m;                 % F=M*A decceleration due to drag in X
AccDragY(1) = -FdragY(1)./ m;                 % F=M*A decceleration due to drag in Y
AccDrag(1,:) = [AccDragX(1) AccDragY(1) 0] 

xImpact_no_Drag = (sind(2*Alpha)*(v0^2))/g
t(5)
for i=2:length(t)
    p(i) = p(i-1) + Vel.*t(i-1) + AccDrag.*(t(i-1).^2);
    
    Vel = Vel + AccDrag.*t(i-1)
    v= sqrt(sum((Vel.^2),"all"))
    
    Velx = v*cosd(Alpha)
    Vely = v*sind(Alpha)- g.*t(i)
    Vel = [Velx Vely 0];        
    
    FdragX = v.*Velx*c;                  % Drag force in X direction
    FdragY = v.*Vely*c;                  % Drag force in Y direction
    AccDragX = -FdragX./ m;                 % F=M*A decceleration due to drag in X
    AccDragY = -FdragY./ m;                 % F=M*A decceleration due to drag in Y
    AccDrag = [AccDragX AccDragY 0]
    if p<0
        break
    end
end
plot(t,p)

Answer 1

[James Tursa]

This makes AccDrag a vector:

AccDrag(1,:) = [AccDragX(1) AccDragY(1) 0]

So here you have a vector on the right hand side and a scalar on the left hand side:

p(i) = p(i-1) + Vel.*t(i-1) + AccDrag.*(t(i-1).^2);

Looks like you could fix this with:

p(i,:) = p(i-1,:) + etc.

BUT, the above equation is wrong. The right hand side needs to be current position + delta position. But you have used formulas that include the absolute time, not the delta time. This will generate an incorrect answer. You need to change how you are updating the position and velocity to use delta times, not absolute times.

Also this test:

if p<0

should be this (if you use the fix above):

if p(i,2)<0