brysonMaxrangeDae(sol,iphase);

function dae = brysonMaxrangeDae(sol,iphase);

global constants

t = sol{1};
x = sol{2};
u = sol{3};
p = sol{4};

xdot = x(:,3).*u(:,1);
ydot = x(:,3).*u(:,2);
vdot = constants.acc-constants.gravity.*u(:,2);
path = u(:,1).^2+u(:,2).^2;
dae = [xdot ydot vdot path];