@vv 

I had copied the lines from my code which contains all the alias statements at the start. Probably should've included that in the code over here, but I thought this conveyed the idea. Could you please elaborate on why alias should be avoided?

Also, I corrected the zeta to Zeta and used this code

alias(beta = beta(t));
alias(Zeta = Zeta(t));
alias(theta = theta(t));

dd := [indets(Ex, specfunc(diff))[]]:
DD := [seq(ddd[i], i = 1 .. nops(dd))]:
EX := eval(Ex, `~`[`=`](dd, DD)):
F := proc (u) options operator, arrow; `if`(limit((eval(u, [beta = Z, Zeta = Z, theta = Z]))*(1/Z), Z = 0) = 0, 0, u) end proc:
Result := eval(map(F, EX), `~`[`=`](DD, dd));

This code returns the correct output. However, when I run it in my worksheet, I get something else.

Ex-Result;
                                                                2
                          3                        3 / d       \ 
-36.86400000 cos(delta__3)  Zeta beta sin(delta__3)  |--- theta| 
                                                     \ dt      / 

   - 16342.88754 sin(psi) beta sin(delta__3) p Zeta

                                                     4   / d  
   + 9.216000000 sin(psi) sin(delta__3) cos(delta__3)  q |--- 
                                                         \ dt 

       \                                             / d       \
  theta| - 5809.338000 cos(psi) beta sin(delta__3) p |--- theta|
       /                                             \ dt      /

                                               / d       \
   + 5805.882000 cos(psi) beta cos(delta__3) q |--- theta|
                                               \ dt      /

                                                / d       \
   + 1.728000000 cos(psi) sin(delta__3) q theta |--- theta|
                                                \ dt      /

                 /  2      \
                 | d       |
   + 2428.922250 |---- beta|
                 |   2     |
                 \ dt      /

                                                         2
   - 3416.512500 V sin(psi) beta sin(delta__3) Zeta theta 

                              3                
   + 16342.88754 cos(psi) beta  sin(delta__3) q

@vv Apologies for that, it's my first time posting here.

The following is a portion of my equation:

Ex := 2428.922250*(diff(beta, t, t))+33032*beta-18.43200000*cos(psi)*beta*(sin(delta__3)^3)*p*r-36.86400000*(cos(delta__3)^3)*Zeta*beta*(sin(delta__3)^3)*((diff(theta, t))^2)-429.7500000*U*cos(psi)*sin(delta__3)*Zeta*Zeta(1, t)+16342.88754*cos(psi)*(beta^3)*sin(delta__3)*q-3416.512500*V*sin(psi)*beta*sin(delta__3)*Zeta*(theta^2)-429.7500000*V*sin(psi)*sin(delta__3)*Zeta*Zeta(1, t)-16342.88754*sin(psi)*beta*sin(delta__3)*p*Zeta+9.216000000*sin(psi)*sin(delta__3)*(cos(delta__3)^4)*q*(diff(theta, t))-5809.338000*cos(psi)*beta*sin(delta__3)*p*(diff(theta, t))-3.456000000*cos(psi)*beta*sin(delta__3)*p*r+149.6160000*cos(psi)*sin(delta__3)*cos(delta__3)*p*Zeta(1, t)+149.6160000*cos(psi)*sin(delta__3)*cos(delta__3)*p*r+5805.882000*cos(psi)*beta*cos(delta__3)*q*(diff(theta, t))+2704.842000*cos(psi)*sin(delta__3)*q*Zeta*r+1.728000000*cos(psi)*sin(delta__3)*q*theta*(diff(theta, t));

I tried using your code, but it returns the same result:

dd := [indets(Ex, specfunc(diff))[]]:
DD := [seq(ddd[i], i = 1 .. nops(dd))]:
EX := eval(Ex, `~`[`=`](dd, DD)):
F := proc (u) options operator, arrow; `if`(limit((eval(u, [beta(t) = Z, zeta(t) = Z, theta(t) = Z]))*(1/Z), Z = 0) = 0, 0, u) end proc;
Result := eval(map(F, EX), `~`[`=`](DD, dd));

I indeed want to retain the derivatives as they are. The assumption that I want to use is that beta, zeta and theta are small, so any term that contains a square or higher power of beta, zeta or theta, or any term that contains a combination of beta, zeta or theta should be neglected. So in this example, the output should be:

Ex := 2428.922250*(diff(beta, t, t))+33032*beta-18.43200000*cos(psi)*beta*(sin(delta__3)^3)*p*r-429.7500000*U*cos(psi)*sin(delta__3)*Zeta*Zeta(1,t)-429.7500000*V*sin(psi)*sin(delta__3)*Zeta*Zeta(1, t)+9.216000000*sin(psi)*sin(delta__3)*(cos(delta__3)^4)*q*(diff(theta, t))-5809.338000*cos(psi)*beta*sin(delta__3)*p*(diff(theta, t))-3.456000000*cos(psi)*beta*sin(delta__3)*p*r+149.6160000*cos(psi)*sin(delta__3)*cos(delta__3)*p*Zeta(1, t)+149.6160000*cos(psi)*sin(delta__3)*cos(delta__3)*p*r+5805.882000*cos(psi)*beta*cos(delta__3)*q*(diff(theta, t))+2704.842000*cos(psi)*sin(delta__3)*q*Zeta*r+1.728000000*cos(psi)*sin(delta__3)*q*theta*(diff(theta, t));

Hope I was able to present it more clearly this time. Thanks for your response!

@vv Added an example :)

The equation looks like the example shown below. Here everything apart from delta__3 is a function of time. Beta, zeta and theta are the main variables (also functions of time), to which I need to apply the small angle approximation. So I need to neglect the terms which contain products of beta, zeta and theta in any combination (ex. beta*theta, zeta⁴ * theta*beta, theta*zeta³ etc) or higher  powers of beta, zeta or theta (ex. beta², zeta⁶, theta⁵ etc).