d = 10;
thetas = find_anglec(d,opt);

figure(1);
clf;
hold on;
for j = 1:2
theta = thetas(j);
v0 = s*[cos(theta); sin(theta)];
tfinal = 2*v0(2)/g;
tref = linspace(0,tfinal)';
y0 = [0; 0; v0];
odeopt = odeset('Events', @hitground);
[tout,yout] = ode45(@(t,y) fball(t,y,opt), tref, y0, odeopt);
plot(yout(:,1), yout(:,2), 'b-');
fprintf('Theta = %g: %e targeting error\n', theta, yout(end,1)-d);
end