function Strela
global g D;
g = 9.81; m=0.2; 	% tíhové zrychlení a hmotnost støely
x0=0; y0=0; v0=30;		% poèáteèní podmínky výstøelu
alfa = pi/6;			% úhel výstøelu
T0 = 0; T1 = 4; 			% èasový interval
cla; hold on;
ind=0; yinit = [x0,v0*cos(alfa),v0*sin(alfa),y0];
for k = 0:0.002:0.004
ind=ind+1;
  C = k; D=C/m;          % souèinitel odporu prostøedí D
  [t,y] = ode45(@f,[T0,T1],yinit);
  plot(y(:,1),y(:,3),'k');
  maxy(ind)=max(y(:,3));
end;
ylim([0 1.05*max(maxy)])
hold off

function z = f(t,y)
  global g D;
  z1 = y(2);
  z2 = -D*y(2)*(y(2)^2+y(4)^2)^(0.5);
  z3 = y(4);
  z4 = -D*y(4)*(y(2)^2+y(4)^2)^(0.5)-g;
  z = [z1;z2;z3;z4];











