Int(sin(a*x)/x,x=0..infinity);
diff(%, a);

I bet that the (relative) small differences are due to numerical rounding errors (it is unclear, what you are using as computational precision). And some changing ordrer of execution (it is not guaranteed in M that it is always the same)

For uploading your sheet: at the end of the menu of the editor (for contributing to this forum) you find a green arrow. Click it. Another menu opens (like in many forums). Use the first step to locate your file on your local machine, then use the upload command, as it is presented. Now choose "insert as link" (not otherwise).

# alias(beta = RootOf(x^3 + x + 1 = 0));
  alias(beta = RootOf(x^3 + x + 1 ,x));  # I prefer that

  expand(q); evala(%);

  expand(t+s); evala(%);

Then both lines work, giving -beta^2-beta;

So you mean the sign of sin(x), x in the Reals - a "square wave"?

de:=rhs(eq1pp);
DE:=convert(de, D);
coeff(DE, D(theta)(t));
#coeff(DE, (D@@2)(theta)(t));

                      Ctheta sin(theta[0] + phi)
                      --------------------------
                              a sin(phi)

Be aware, that 'coefficients' - like for polynomials - may
need to expand the expression first.

PS: do not hesitate to use line breaks for longer formulae :-)

convert(%, expln) gives an answer as well.

Using 'gfun' it turns out to be a root of 56*x^3 - 42*x - 13

'proof':

  RootOf(56*x^3-42*x-13, index=1) / cos(arccos(13/14)/3) = 1;
  allvalues(%); convert(%, expln); evala(%); is(%);

                    true

g:= x -> BesselJ(0,x);
plot(g(x), x = -4 .. 4);

Summarising the discussion and argumentations in a sheet ...

MP_int_products_sinc.mws
MP_int_products_sinc.pdf

DTEdt:=aE*qE*A/(alpha+TR_SS)-muE*TE+rE TE   (1-TE/TEmax);

DTEdt:=aE*qE*A/(alpha+TR_SS)-muE*TE+rE*TE * (1-TE/TEmax);

You may want to use 'else' as last statement. I even type it always before writing 'end if' in any case - before filling the body, just as a frame.

For clarity I prefer a 'return' in the if-branches (not needed), it depends on the 'language' you use and not all programming languages behave the same way on constructions.

A suggestion for clarity: do not use the same name beta for different things.

If you mean a function then name it p, for example, like p(0.787*beta)

Then it makes it clear, that you have to tell us what p stands for,
it does not make much sense to say 'can be any value'

Those sheets in Standard Mode are very unhandily, so I use the better interface :-)

Find a sheet attached, the value for i=53 is 0.606852403708179e-2

MP_2dim-Int-numerica.mws

One may try to investigate closer for x0:=12236889.2928742 (having Digits=15) due to

  eps:=1e-5;
  m/10^24;
  plot([Re(%), Im(%)], x= x0*(1-eps) .. x0*(1+eps), color=[red,blue]);

I suggest that, after studying asympt(m, x, 1) and piecwise and some scaling, 
v:= m/( (1000 + 10000*I)*x^3*exp(0.4e-7*x) ) and some 'Brent Zero Finder'

Try a polynomial, degree = 7 (since there are 8 equations), now solve it

If have a Standard sheet or you open as such, then you can export as pdf.

You almost always can save as html, which can be handled by Word, where
input can be maintained, outputs a graphics.

And you generally can install some pdf or ps printer on your machine, which
can print as device anyway.