@Kaare112 We may assume that F_D(v) is supposed to return the air resistance force on input of the velocity (here a list). If so, there is a missing minus sign since that force is in the opposite direction of the velocity.
Thus it should be
k:=-evalf((3*Pi*mu*d+(Pi/8)*C_0*rho*(d^2)*sqrt(vx^2+vy^2))/m); #notice minus added
Also it appears that Eskritt is one Euler step in solving dv/dt = F/m.
Thus to do several steps in solving that ode you could do:
#With the first two lines from your code above:
phi:=Pi/5: v0:=evalf([v_0*cos(phi), v_0*sin(phi)]); v[0]:=v0;
h:=0.002: N:=150:
#and then
V:=Vector(N); #V will be a vector having velocity (lists) as entries
v_0:=10: #Initial speed set arbitrarily here for illustration
V[1]:=v0:
for i from 2 to N do V[i]:=Eskritt(V[i-1],h) end do:
V[1..10];
V[21..30];
V[N-9..N];
If you want to see the whole vector do
interface(rtablesize=infinity);
V;
plot(V,labels=["vx","vy"]);
Now that you have the velocity you should be able to find the position.
But maybe both should be found at the same time, i.e. corresponding to the system of 4 odes
dv/dt = F/m
dr/dt = v
both vector equations.

Since you have already done some work you will have more luck getting useful responses if you show us what you have done.What are the assumptions: No air resistance?
Do you have to use Euler's method?

So your effort so far consists in posting an image in the questions section of MaplePrimes?

Actually you can also use Heaviside via piecewise like this:

a0(t):=piecewise(x(t)+a[1](t)>0,-x(t),a[1](t));
a0(t):=convert(a0(t),Heaviside);
This has the same effect of making it possible for convertsys to do its job of isolating the highest derivatives.

Actually you can also use Heaviside via piecewise like this:

a0(t):=piecewise(x(t)+a[1](t)>0,-x(t),a[1](t));
a0(t):=convert(a0(t),Heaviside);
This has the same effect of making it possible for convertsys to do its job of isolating the highest derivatives.

Since you found no problem in Maple with the parameter values mentioned maybe you could give us an example of a pair of parameter values zeta and n for which problems do exist in Maple.

@acdah The way you proceed it does not make sense to assign values to theta, r, epsilon.
That was  really meant by saying you were stuck.

If you want approximations they should take place before integration. So approximate the integrand.

@acdah  Being of type numeric means being actual numbers, i.e. integers, rationals, or floating point numbers. So you cannot change the types of the variables, but you can assign values of type numeric to them, as in alpha:=1.2345;

You have an unevaluated integral in mH from the integral in mean_integral, where you have something like
Int( f(r,alpha,beta,epsilon,theta,M), theta=0..Pi,method= _Gquad).
_Gquad is a numerical method and would require r,alpha,beta,epsilon,M all to be of type numeric, but in your expression they are names.
If you replace the definition by
mean_integral := 2*Pi*int(simplify(mean_square)*area_element, theta = 0 .. Pi);
you will see that it returns unevaluated. I think you are stuck there.

@N00bstyle Your procedure reductiefactorvoetgangershoogte is unnecessarily complicated and involves a lot of calls to eigenfrequentieigvvoetgangers. But in this case the main problem is that evalb (which is used in the if statements) cannot handle Pi:
try
if Pi<3 then 7 else 0 end if;
evalb(Pi<3);
You could use 'is' instead of evalb as in
if is(Pi<3) then 7 else 0 end if;
is(Pi<3);
Instead you could replace Pi by evalf(Pi).
So a suggestion for your procedure reductiefactorvoetgangershoogte would be to change the procedure to
reductiefactorvoetgangershoogte := proc (_hoogte)
local eg; eg := evalf(eigenfrequentieigvvoetgangers(_hoogte));
if eg <= 1.25 then 0 elif 1.25 < eg and eg < 1.7 then (eg-1.25)/(1.7-1.25) elif 1.7 <= eg and eg <= 2.1 then 1 elif 2.1 < eg and eg < 2.3 then (2.3-eg)/(2.3-2.1) elif 2.3 <= eg and eg <= 2.5 then 0 elif 2.5 < eg and eg < 3.4 then .25*(eg-2.5)/(3.4-2.5) elif 3.4 <= eg and eg <= 4.2 then .25 elif 4.2 < seg and seg < 4.6 then .25*(4.6-eg)/(4.6-4.2) elif 4.6 <= eg then 0 end if end proc;
(I'll let you do the indenting)
My points are: eigenfrequentieigvvoetgangers(_hoogte) is evaluated only once and then evalf is aplied.
However, I still think you should use unapply and piecewise instead of your present approach.
Why insist on making it complicated?

@N00bstyle Your procedure reductiefactorvoetgangershoogte is unnecessarily complicated and involves a lot of calls to eigenfrequentieigvvoetgangers. But in this case the main problem is that evalb (which is used in the if statements) cannot handle Pi:
try
if Pi<3 then 7 else 0 end if;
evalb(Pi<3);
You could use 'is' instead of evalb as in
if is(Pi<3) then 7 else 0 end if;
is(Pi<3);
Instead you could replace Pi by evalf(Pi).
So a suggestion for your procedure reductiefactorvoetgangershoogte would be to change the procedure to
reductiefactorvoetgangershoogte := proc (_hoogte)
local eg; eg := evalf(eigenfrequentieigvvoetgangers(_hoogte));
if eg <= 1.25 then 0 elif 1.25 < eg and eg < 1.7 then (eg-1.25)/(1.7-1.25) elif 1.7 <= eg and eg <= 2.1 then 1 elif 2.1 < eg and eg < 2.3 then (2.3-eg)/(2.3-2.1) elif 2.3 <= eg and eg <= 2.5 then 0 elif 2.5 < eg and eg < 3.4 then .25*(eg-2.5)/(3.4-2.5) elif 3.4 <= eg and eg <= 4.2 then .25 elif 4.2 < seg and seg < 4.6 then .25*(4.6-eg)/(4.6-4.2) elif 4.6 <= eg then 0 end if end proc;
(I'll let you do the indenting)
My points are: eigenfrequentieigvvoetgangers(_hoogte) is evaluated only once and then evalf is aplied.
However, I still think you should use unapply and piecewise instead of your present approach.
Why insist on making it complicated?

@N00bstyle If after having defined eigenfrequentie in your version, you try the command
eigenfrequentie;
then you reveive the error
Error, (in qigvvoetgangers) cannot determine if this expression is true or false: 0 <= 0.153846153846154e-1*(65.0*hoogte+.967200000000000)/hoogte-21/20
That is where you have to ask yourself what to do. In fact you must have noticed the problem yourself (I guess) since you don't use eigenfrequentie in the plot command defining ploteigenfrequentie, but rather the expression used to define eigenfrequentie.

@N00bstyle If after having defined eigenfrequentie in your version, you try the command
eigenfrequentie;
then you reveive the error
Error, (in qigvvoetgangers) cannot determine if this expression is true or false: 0 <= 0.153846153846154e-1*(65.0*hoogte+.967200000000000)/hoogte-21/20
That is where you have to ask yourself what to do. In fact you must have noticed the problem yourself (I guess) since you don't use eigenfrequentie in the plot command defining ploteigenfrequentie, but rather the expression used to define eigenfrequentie.

@Markiyan Hirnyk Try

inttrans[invlaplace](s^(-2.3)/(s+c),s,t); #No problem
inttrans[invlaplace](s^(-2.3)*cos(s)/(s+c),s,t); #Not found
That the last one is not found is just because invlaplace isn't perfect.

In the question p was assumed constant. However, it may be worth mentioning that the non-constant case can be handled as well, but we must use DEtools[mult]:

with(DEtools):
L:=(Dx-p(x))^2;
diffop2de(L,f(x),[Dx,x]); # does not correspond to composition
#This gives the correct version:
mult(Dx-p(x),Dx-p(x),[Dx,x]);
diffop2de(%,f(x),[Dx,x]);