my_velocity.m

function [vx,vy,vz] = my_velocity(x,y,z,t)
    % Help documentation of "my_velocity"
    % This function computes the velocity in the x- and y-directions given
    % three vectors of position in x- and y-directions as a function of time
    % x = x-position
    % y = y-position
    % z = z-position
    % t = time
    % output
    % vx = velocity in x-direction
    % vy = velocity in y-direction
    % vz = velocity in z-direction

    vx=zeros(length(t),1);
    vy=zeros(length(t),1);
    vz=zeros(length(t),1);
    
    vx(1:end-1) = diff(x)./diff(t); % calculate vx as delta x/delta t
    vy(1:end-1) = diff(y)./diff(t); % calculate vy as delta y/delta t
    vz(1:end-1) = diff(z)./diff(t); % calculate vy as delta y/delta t
    
    vx(end) = vx(end-1);
    vy(end) = vy(end-1);
    vz(end) = vz(end-1);
    
end