I consider it inappropriate to start your question with "In the name of God". Who's god? The Christian god or what? It doesn't matter. What has your question got to do with God?

You are sure to antagonize Christians if they think you might be a Muslim, Muslims if they think you are Christian, and atheists in any case. Thus you are less likely to get responses.

You are more likely to get good answers if you provide more details, in this case the integral in question.

Here are two animates inside each other. I didn't try three.

animate(animate,[implicitplot3d,[eval(f,{a=1/2})=50, x=0..40 ,  y=0..40, z=0..40,axes=boxed,caption=typeset("c = ",c,", b= ",b)],b=0..1,frames=7,paraminfo=false],c=0..1,frames=7,paraminfo=false);

Here are two animates inside each other. I didn't try three.

animate(animate,[implicitplot3d,[eval(f,{a=1/2})=50, x=0..40 ,  y=0..40, z=0..40,axes=boxed,caption=typeset("c = ",c,", b= ",b)],b=0..1,frames=7,paraminfo=false],c=0..1,frames=7,paraminfo=false);

I'm guessing that this is a physics problem of some sort, KE being kinetic energy and PE being potential energy.

Is PI meant to be Pi (i.e. our familiar constant Pi)? Not that it helps though. The integral seems difficult.

In the situation could the integrals be definite integrals, then at least numerical computation could be done for specific values of the constants?

I'm guessing that this is a physics problem of some sort, KE being kinetic energy and PE being potential energy.

Is PI meant to be Pi (i.e. our familiar constant Pi)? Not that it helps though. The integral seems difficult.

In the situation could the integrals be definite integrals, then at least numerical computation could be done for specific values of the constants?

You could use ghostscript.

http://pages.cs.wisc.edu/~ghost/

You could use ghostscript.

http://pages.cs.wisc.edu/~ghost/

OK, something like the result of this?

int(int(exp(x^3)*sin(x)*y,y),x);

Your description of the problem is not simple enough for me to understand. I suggest that you post the code.

You may have noticed that nothing is plotted in the left half plane.

x^(1/3) is the principal value of the complex cube root. To get the one intended use surd(x,3) instead:

dfieldplot(diff(y(x), x) = sin(x*y(x))+3*x^2+surd(x,3)+4*y(x)^2+5*y(x), y(x), x = -5 .. 5, y = -5 .. 5, title = `tangent field`, color = black);

You may have noticed that nothing is plotted in the left half plane.

x^(1/3) is the principal value of the complex cube root. To get the one intended use surd(x,3) instead:

dfieldplot(diff(y(x), x) = sin(x*y(x))+3*x^2+surd(x,3)+4*y(x)^2+5*y(x), y(x), x = -5 .. 5, y = -5 .. 5, title = `tangent field`, color = black);

The equation I called E1, i.e.

diff(y(x), x) = sqrt(1-y(x)^2)/sqrt(1+x^2);

does not satisfy the usual requirements for uniqueness of solutions at points (x0, 1) or (x0, -1) for any x0.

Uniqueness of solutions at (x0, y0) is guranteed for the ode y' = f(x, y) with y(x0) = y0, if f and its partial derivative wrt. y are continuous in a neighborhood of (x0, y0).

These requirements are not satisfied for E1 for any point (x0, 1) or (x0, -1) since the partial derivative wrt. y of 

sqrt(1-y^2)/sqrt(1+x^2)

doesn't even exist at any such point.

Now the requirements are sufficient (not necessary), but there are in fact infinitely many solutions passing through any point (x0, 1), and the same holds for (x0,-1).

Your original equation E gives rise to two odes, which I called E1 and E2. Solutions to E1 has nondecreasing solutions and solutions to E2 have nonincreasing solutions. This allows for infinitely many oscillating solutions all starting with y(0) = 0.

The equation I called E1, i.e.

diff(y(x), x) = sqrt(1-y(x)^2)/sqrt(1+x^2);

does not satisfy the usual requirements for uniqueness of solutions at points (x0, 1) or (x0, -1) for any x0.

Uniqueness of solutions at (x0, y0) is guranteed for the ode y' = f(x, y) with y(x0) = y0, if f and its partial derivative wrt. y are continuous in a neighborhood of (x0, y0).

These requirements are not satisfied for E1 for any point (x0, 1) or (x0, -1) since the partial derivative wrt. y of 

sqrt(1-y^2)/sqrt(1+x^2)

doesn't even exist at any such point.

Now the requirements are sufficient (not necessary), but there are in fact infinitely many solutions passing through any point (x0, 1), and the same holds for (x0,-1).

Your original equation E gives rise to two odes, which I called E1 and E2. Solutions to E1 has nondecreasing solutions and solutions to E2 have nonincreasing solutions. This allows for infinitely many oscillating solutions all starting with y(0) = 0.

@BSchor You may try

sin(N*Pi) assuming N::integer;

and also sin(N*Pi);

The procedure sin knows how to handle Pi, try

showstat(sin);

I don't believe simplify is involved.

Another point, try

sum(k^2,k=1..N);
eval(%,N=3/2);
add(k^2,k=1..3/2);
sum(k^2,k=1..3/2);

You see that sum tries to find a general formula in the first case. That formula makes sense even for values of N that are not integers.