How to define the function for 'velocityspinsolver' in ode45 refer to the scipt below

I attached this with the journal. i know the 'velocityspinsolver' is refer to this equations:

Here the fully original script, but i cannot run because i don't know to write that 'velocityspinsolver':

%Program intro disp
disp('This program was written to calculate the fastest volleyball serve');
disp('given several user-defined values. The program assumes the weight');
disp('of the ball, the height of the net, the size of the ball, and the');
disp('drag coefficient of the the ball. Answer questions in feet');

%Input parameters
height = input('What height is the volleyball served from: \n'); 
phi= input('What angle are you serving with respect to the sideline (in degrees): \n'); 
dtonet2 = input('What is the shortest distance to the net from where you are standing: \n'); 

dfromnet2 = input('How far from the net would you like it to land: \n'); 
angularvelocity2 = input('How much spin is put on the ball (in revolutions per second): \n');

%Define assumed parameters
netheight = 7.6919; %in ft, includes ball size 
global dragcoef; 
dragcoef = 0.35; 
global densityofair; 
densityofair = 0.0729; %temp =70, altitude = 600 ft, in lb/ftˆ3 
global areaball; 
areaball = 0.3808; %in sqare ft 
global massball;
massball = 0.59375; %in pounds 
global gravity; 
gravity = -32.174;%ft/s/s 
global radius;%ft 
radius = 0.4; 
global angularvelocity; %???? 
angularvelocity=angularvelocity2*2*pi; 
global spincoef; 
spincoef = 0.3;

%Calculate true distance from and to net
phiradians = phi*pi/180; 
dtonet = dtonet2/cos(phiradians); 
dfromnet = dfromnet2/cos(phiradians);

%Initial loop, moves large incriments
q=0; 
for theta=0:(pi/90):(pi/2) 
    for speed=150:-5:0 
        V0=[speed*cos(theta) speed*sin(theta)];
        [time,velocities] = ode45('velocityspinsolver',0:.005:5, V0); 
        X = zeros (length(time),2); 
        X(1,2)=height; 
        for i = 1:(length(time)-1) 
            X(i+1,1)=(X(i,1)+(time(i+1)-time(i))*velocities(i,1)); 
            X(i+1,2)=(X(i,2)+(time(i+1)-time(i))*velocities(i,2));
        end
        findnet=(X(:,1)<dtonet); 
        if sum(findnet)>length(time)
            if X(sum(findnet)+1,2)> netheight 
                b=sum(X(:,2)>0); 
                if (X(b,1)>dtonet) && (X(b,1)<(dtonet+dfromnet)) 
                    q=1; 
                end
            end
        end
        if q==1 
            break
        end
    end
    if q==1 
       break
    end
end

%Fine tuned loop: This loop uses the results from the initial loop to
%obtain a result accurate to within a degree and a foot per second
q=0; 
for theta2=(theta-(pi/10)):(pi/180):(pi/2)
for speed2=150:-1:(speed-20)
V0=[speed2*cos(theta2) speed2*sin(theta2)];
[time,velocities] = ode45('velocityspinsolver',0:.005:5, V0);
X = zeros (length(time),2);
X(1,2)=height;
for i = 1:(length(time)-1)
X(i+1,1)=(X(i,1)+(time(i+1)-time(i))*velocities(i,1));
X(i+1,2)=(X(i,2)+(time(i+1)-time(i))*velocities(i,2));
end
findnet=(X(:,1)<dtonet);
if sum(findnet)> length(time)
if X(sum(findnet)+1,2)> netheight
b=sum(X(:,2)>0);
if (X(b,1)>dtonet) && (X(b,1)<(dtonet+dfromnet))
q=1;
end
end
end
if q==1
break
end
end
if q==1
end
end
%Displays results
disp('Ball should be hit at a speed of') 
disp(speed2) 
disp('and an angle of ') 
disp(theta2*180/pi) 
disp(' degrees.') 
disp('A plot of the trajectory of the ball is displayed.') 
plot(X(:,1),X(:,2))

[Lance Anil]

Same issue. Were you able to solve it?