function out=waterwheel_nd(t,all_vector,p)

% extract ODE parameter values from input vector
m         = all_vector(1:p.N);  % grab mass vector from input
theta1    = all_vector(p.N+1);  % theta for cup #1
theta     = p.thetadist+theta1; % map theta1 onto all cups
omega     = all_vector(end);    % grab omega vector from input

% nondimensional mass as a function of nondimensional time (diff eq)         
dmdt=real(-p.D1*sqrt(m)+p.D2*sin(theta/2).^p.large_n.*(m<1));

% nondimensional theta as a function of nondimensional time (diff eq)
dtheta1dt = omega;

% nondimensional omega as a function of nondimensional time (diff eq)
domegadt=(-sum(m.*sin(theta))-p.D3*omega-omega*sum(dmdt))/sum(m);

% combine output for ode output
out       = [dmdt; dtheta1dt; domegadt];
