Нахождение максимальной высоты от угла выстрела

> restart;
> g := 9.81; mu := 0.182e-1; r := .1; c := .4; rho := 1.2928; m := 4; h0 := 2.5; v0 := 10; a0 := (1/18)*Pi; an := (1/3)*Pi; n := 7; astep := (an-a0)/n; tmax := 5; k1 := (1/6)*Pi*r*mu; k2 := (1/2)*r^2*c*rho*Pi;
> for i from 0 to n do A0[i] := a0+i*astep; dsol := dsolve({diff(x(t), t, t) = -(k1+k2*sqrt((diff(x(t), t))^2+(diff(y(t), t))^2))*(diff(x(t), t))/m, diff(y(t), t, t) = -(k1+k2*sqrt((diff(x(t), t))^2+(diff(y(t), t))^2))*(diff(y(t), t))/m-g, x(0) = 0, y(0) = h0, (D(x))(0) = v0*cos(A0[i]), (D(y))(0) = v0*sin(A0[i])}, numeric, output = listprocedure); X[i] := eval(x(t), dsol); Y[i] := eval(y(t), dsol); Vx[i] := eval(diff(x(t), t), dsol); Vy[i] := eval(diff(y(t), t), dsol); tgrd[i] := fsolve(Y[i](t), t = 0 .. tmax) end do;
> plot({seq([X[i](t), Y[i](t), t = 0 .. tgrd[i]], i = 0 .. n)});
> plot({seq([t, sqrt(Vx[i](t)^2+Vy[i](t)^2), t = 0 .. tgrd[i]], i = 0 .. n)})
> Xa := evalf([seq([A0[i], X[i](tgrd[i])], i = 0 .. n)]); with(CurveFitting); IP := proc (x) options operator, arrow; PolynomialInterpolation(Xa, x) end proc; IP(x)
> plot({Xa, IP(a)}, a = a0 .. an)
> Amax := fsolve(diff(IP(a), a) = 0, a = a0 .. an); Xmax := IP(Amax)