Plotting projectile motion with variable drag

I have an assignment to plot the motion of a projectile as it travels a distance, ideally to (50,30), 'catches' a drone, and then falls to the ground with a parachute. This means the mass, size of the projectile, and coefficient of drag change with distance. The first part of the assignment asks to just plot the motion given an angle using three functions: stateDeriv, stepRungeKutta, and ivpSolver I have attempted to make a final function dragForce to calculate drag at a given point. The initial velocity is 50m/s and it starts at (0,1).

function  [Fx,Fy]  = DragForce(z)
% Function to find the combined drag force at various stages
%   given z = [Xx; Xy; Vx; Vy] find the drag force in x and y
Vx = z(3);      % velocity in x (m/s)
Vy = z(4);      % velocity in y (m/s)
theta = 60;     % starting angle (degrees)
rho = 1.225;    % density of air (kg/m^3)
g = 9.81;       % gravitational force (N)
x = z(1);       % displacement in x (m)
Y = z(2);       % displacement in y (m)

      if x < 50
        M = 0.5;
        A = 0.0005;
        C = 0.1;

    else  
        M = 1;
        A = 1.5;
        C = 0.9;

    end

     Fx = (0.5 * C * A * rho * (Vx)^2)/M ;

     Fy =  (0.5 * C * A * rho * (Vy)^2)/M;
     end

stateDeriv uses coupled first order ODE's to find the state derivative of the projectile

    function [dz] = stateDeriv(x,z,Fx,Fy)
     % Calculate the state derivative for the motion of the projectile

    %     DZ = stateDeriv(X,Z) computes the derivative DZ = [Vx Vy; Ax Ay] of the 
    %     state vector Z = [Xx; Vx; Xy; Vy], where X is displacement, V is velocity,
    %     and A is acceleration.

    g = 9.81; % gravitational force kg/m^3

    [Fx,Fy] = DragForce(z);

% You have to complete the next two lines to calculate the state 
% derivatives (Hint: draw the free-body diagram and use Newton's 2nd Law

    dz1 = z(2);
    dz2 = Fx;

    dz3 = z(4);
    dz4 = Fy - g;

    dz = [dz1; dz3; dz2; dz4];
    end

stepRungeKutta shouldn't be edited (apparently) but I replaced t and dt with x and dx because I want to plot y against x

    function znext = stepRungeKutta(x,z,dx)
    % stepRungeKutta    Compute one step using the RungeKutta method
    % 
    %     ZNEXT = stepRungeKutta(T,Z,DT) computes the state vector ZNEXT at the next
    %     time step T+DT

    % Calculate the state derivative from the current state

    A = dx * stateDeriv(x, z);
    B = dx * stateDeriv(x + (dx/2), z + (A/2));
    C = dx * stateDeriv(x + (dx/2), z + (B/2));
    D = dx * stateDeriv(x + dx, z + C);

    % Calculate the next state vector from the previous one using Runge Kutta 4
    znext = z + (A + (2 * B) + (2 * C) + D) / 6; 
    end

ivpSolver I am completely stuck on. I think that my previous functions are okay, but this is where I have issues.

    function [x,z] = ivpSolver(x0,z0,dx,xend)
    % ivpSolver    Solve an initial value problem (IVP) and plot the result
    % 
    %     [x,Z] = ivpSolver(T0,Z0,DX,XE) computes the IVP solution using a step 
    %     size Dx, beginning at distance x0 and initial state Z0 and ending at distance 
    %     XEND. The solution is output as a time vector T and a matrix of state 
    %     vectors Z.

    % Set initial conditions
    x = x0
    g = -9.81;    % Acceleration due to gravity.
    z(:) = z0 

    % Continue stepping until the end time is exceeded
    n=1;
    while x(n) <= xend
      % Increment the time vector by one time step
      x(n+1) = x(n) + dx;

      % Apply Runge Kutta's method for one time step
      z(:,n+1) = stepRungeKutta(x(n), z(:,n), dx);

      n = n+1;
  end

% Compute acceleration from velocity (not needed for solving the
% problem, but nice to have for completeness)
ddz = diff(z(2,:)) / dx;

% Plot the result
figure(1)
plot(x,z(1,:),'b')

xlim([x0 xend])
xlabel('x distance, m')
ylabel('y distance, m')
legend('Displacement')

These are my four functions, I'm getting these errors: Index exceeds matrix dimensions.

Error in DragForce (line 4) Vx = z(3); % velocity in x (m/s)

Error in stateDeriv (line 9) [Fx,Fy] = DragForce(z);

Error in stepRungeKutta (line 9) A = dx * stateDeriv(x, z);

Error in ivpSolver (line 25) z(:,n+1) = stepRungeKutta(x(n), z(:,n), dx);