@Alina8faq I now see what you did. You were in 2D input mode and then punched sin^2 resulting in

after that you moved the cursor down (with the arrow key) and punched (a+b).
One advantage of using 1D-input (as I always do) is that descriptions like the short one I gave above are totally unnecessary. You just write characters one by one and your commands may as well be written in a text editor. So when you wrote
K:=(sin^2(a+b)-sin^2(a-b))
then you didn't just punch those characters you also moved the cursor when necessary.

If you want you can change to 1D-input (Maple input) and worksheet mode under the menu items
Tools/Options/Display and Tools/Options/Interface, respectively.

@Alina8faq I now see what you did. You were in 2D input mode and then punched sin^2 resulting in

after that you moved the cursor down (with the arrow key) and punched (a+b).
One advantage of using 1D-input (as I always do) is that descriptions like the short one I gave above are totally unnecessary. You just write characters one by one and your commands may as well be written in a text editor. So when you wrote
K:=(sin^2(a+b)-sin^2(a-b))
then you didn't just punch those characters you also moved the cursor when necessary.

If you want you can change to 1D-input (Maple input) and worksheet mode under the menu items
Tools/Options/Display and Tools/Options/Interface, respectively.

@Alina8faq If you use 1D input (aka Maple input), then what you wrote results in 0:
K:=(sin^2(a+b)-sin^2(a-b));
returns 0. The reason is that 2(a+b) is interpreted as the constant function 2 applied to a+b, thus 2(a+b) = 2. In the same way 2(a-b) = 2. Thus (sin^2(a+b)-sin^2(a-b)) evaluates to sin^2 - sin^2 , i.e. 0.
If I insert K:=(sin^2(a+b)-sin^2(a-b)) as 2D input:

and then do
lprint(K);
then I get the following printout
sin^2*(a+b)-sin^2*(a-b)
and that is surely not what you want either, because sin^2*(a+b) is interpreted as the function
sin^2 times (a+b).

@Alina8faq If you use 1D input (aka Maple input), then what you wrote results in 0:
K:=(sin^2(a+b)-sin^2(a-b));
returns 0. The reason is that 2(a+b) is interpreted as the constant function 2 applied to a+b, thus 2(a+b) = 2. In the same way 2(a-b) = 2. Thus (sin^2(a+b)-sin^2(a-b)) evaluates to sin^2 - sin^2 , i.e. 0.
If I insert K:=(sin^2(a+b)-sin^2(a-b)) as 2D input:

and then do
lprint(K);
then I get the following printout
sin^2*(a+b)-sin^2*(a-b)
and that is surely not what you want either, because sin^2*(a+b) is interpreted as the function
sin^2 times (a+b).

@fairuzalwani I tried the following, where for my two examples I need simplify(z-ic) in order for the loop to stop. In your case it will depend on the kind of map (Map) you got.
I added a third example below.

Orbit:=proc(Map,ic,Nmax::posint:=100)
  local orbit,z,t;
  orbit:=Vector();
  orbit(1):=ic;
  z:=ic;
  for t from 2 to Nmax while simplify(z-ic)<>0 or t=2 do
    z:=Map(z);
    orbit(t):=z
  end do;
  orbit
end proc;

##########################################
f:=z->z^3;
Orbit(f,exp(I*2*Pi/5),10);
orb:=simplify(%);
plots:-complexplot(exp(I*t),t=0..2*Pi,color=blue):
plots:-complexplot(convert(orb,list),style=point,symbolsize=20,symbol=solidcircle):
plots:-display(%%,%);
#Try the same with only the word 'simplify' removed in the procedure Orbit
##########################################
f:=z->exp(I*Pi/10)*z;
Orbit(f,exp(I*2*Pi/5));
orb:=simplify(%);
plots:-complexplot(exp(I*t),t=0..2*Pi,color=blue):
plots:-complexplot(convert(orb,list),style=point,symbolsize=20,symbol=solidcircle):
plots:-display(%%,%);
#Again you may try the same with only the word 'simplify' removed in the procedure Orbit
####################
#Third example:
If your function maps a set of lists of pairs of integers into itself then the simple check z<>ic should do:
Orbit:=proc(Map,ic,Nmax::posint:=100)
  local orbit,z,t;
  orbit:=Vector();
  orbit(1):=ic;
  z:=ic;
  for t from 2 to Nmax while z<>ic or t=2 do
    z:=Map(z);
    orbit(t):=z
  end do;
  orbit
end proc;
##Simple example
f:=L->modp~(L+[8,9],5);
f([1,1]);
Orbit(f,[1,1],10);

@fairuzalwani I tried the following, where for my two examples I need simplify(z-ic) in order for the loop to stop. In your case it will depend on the kind of map (Map) you got.
I added a third example below.

Orbit:=proc(Map,ic,Nmax::posint:=100)
  local orbit,z,t;
  orbit:=Vector();
  orbit(1):=ic;
  z:=ic;
  for t from 2 to Nmax while simplify(z-ic)<>0 or t=2 do
    z:=Map(z);
    orbit(t):=z
  end do;
  orbit
end proc;

##########################################
f:=z->z^3;
Orbit(f,exp(I*2*Pi/5),10);
orb:=simplify(%);
plots:-complexplot(exp(I*t),t=0..2*Pi,color=blue):
plots:-complexplot(convert(orb,list),style=point,symbolsize=20,symbol=solidcircle):
plots:-display(%%,%);
#Try the same with only the word 'simplify' removed in the procedure Orbit
##########################################
f:=z->exp(I*Pi/10)*z;
Orbit(f,exp(I*2*Pi/5));
orb:=simplify(%);
plots:-complexplot(exp(I*t),t=0..2*Pi,color=blue):
plots:-complexplot(convert(orb,list),style=point,symbolsize=20,symbol=solidcircle):
plots:-display(%%,%);
#Again you may try the same with only the word 'simplify' removed in the procedure Orbit
####################
#Third example:
If your function maps a set of lists of pairs of integers into itself then the simple check z<>ic should do:
Orbit:=proc(Map,ic,Nmax::posint:=100)
  local orbit,z,t;
  orbit:=Vector();
  orbit(1):=ic;
  z:=ic;
  for t from 2 to Nmax while z<>ic or t=2 do
    z:=Map(z);
    orbit(t):=z
  end do;
  orbit
end proc;
##Simple example
f:=L->modp~(L+[8,9],5);
f([1,1]);
Orbit(f,[1,1],10);

@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.