theta = pi/3;
v0 = s*[cos(theta); sin(theta)];

% Compute the reference trajectory (absent air resistance)
% x(y) = s*cos(theta)*t
% y(t) = -g*t^2/2 + s*sin(theta)*t
% up to time tfinal = 2*s*sin(theta)/g
%
tfinal = 2*v0(2)/g;
tref = linspace(0,tfinal)';
xref = v0(1)*tref;
yref = (v0(2)-g/2*tref).*tref;

% Compute the same reference trajectory with ode45 (with drag)
y0 = [0; 0; v0];
odeopt = odeset('Events', @hitground);
[tout,yout] = ode45(@(t,y) fball(t,y,opt), tref, y0, odeopt);

% Do a comparison visually
plot(xref, yref, 'r:', yout(:,1), yout(:,2), 'b-');