@nMaple  Here is an animation with front moving part filled with yellow color (the initial yellow region bounded by y=x^2+1, x=0, and y=5 and rotation about x=-1):

restart;

with(plots): with(plottools):

ax := spacecurve([-1, t, 0], t = -1 .. 6, color = red, thickness = 2, linestyle = dash):

s := plot3d([x, y, 0], x = 0 .. 2, y = x^2+1 .. 5, style = surface, color = yellow):

F := q->rotate(q, alpha, [[-1, 0, 0], [-1, 1, 0]]):

s0 := animate(display, [F(s)], alpha = 0 .. 2*Pi, frames = 60):

s1 := animate(plot3d, [[(x+1)*cos(t)-1, x^2+1, -(x+1)*sin(t)], x = 0 .. 2, t = 0 .. alpha, style = surface, color = blue], alpha = 0 .. 2*Pi, frames = 60):

s2 := animate(plot3d, [[(x+1)*cos(t)-1, 5, -(x+1)*sin(t)], x = 0 .. 2, t = 0 .. alpha, style = surface, color = blue], alpha = 0 .. 2*Pi, frames = 60):

s3 := animate(plot3d, [[cos(t)-1, y, -sin(t)], y = 1 .. 5, t = 0 .. alpha, style = surface, color = blue], alpha = 0 .. 2*Pi, frames = 60):

display([s0, ax, seq(s || i, i = 1 .. 3), s], scaling = constrained, lightmodel = light3, axes = normal, orientation = [60, 70], orientation = [70, 75], labels = [x, y, z]);

@nMaple 

restart;

with(plots):  with(plottools):

ax := spacecurve([-1, t, 0], t = -1 .. 6, color = red, thickness = 2, linestyle = dash):

s1 := animate(plot3d, [[(x+1)*cos(t)-1, x^2+1, -(x+1)*sin(t)], x = 0 .. 2, t = 0 .. alpha, style = surface, color = blue], alpha = 0 .. 2*Pi, frames = 100):

s2 := animate(plot3d, [[(x+1)*cos(t)-1, 5, -(x+1)*sin(t)], x = 0 .. 2, t = 0 .. alpha, style = surface, color = blue], alpha = 0 .. 2*Pi, frames = 100):

s3 := animate(plot3d, [[cos(t)-1, y, -sin(t)], y = 1 .. 5, t = 0 .. alpha, style = surface, color = blue], alpha = 0 .. 2*Pi, frames = 100):

display(ax, seq(s||i, i = 1 .. 3), scaling = constrained, lightmodel = light3, axes = normal, orientation = [60, 70], labels = [x, y, z]);

See my answer in this thread

@Markiyan Hirnyk  I didn't see your comment before posting mine. In any case, a visualization is always helpful as it gives confidence that there are no other solutions.

@patient   of course. It is not necessary to write  method=rkf45, since this method is a method by default. 

restart;

Eq:=diff(u1(z),z)=9/4*(84*(-z)^(11/2)-8*z^7+540*(-z)^(5/2)+324*z^4-324*z)*u1(z)/((z^2+3*sqrt(-z))*(2*(-z)^(3/2)+3)*(8*z^6+66*(-z)^(9/2)-189*z^3+216*(-z)^(3/2)+81)):

Sol:=dsolve({Eq, u1(-2)=1}, numeric, output=listprocedure);

u1:=rhs(Sol[2]);

u0(z):=z^2/3+sqrt(-z):

U:=unapply(eval(u0(z)+t*u1(z), z=x/t),x,t);

plot(U(x,1), x=-2..0, color=red, thickness=2); 

@Markiyan Hirnyk  

In your example, a significant difference between the volume and the number of integer points can be explained by the small  radius of the ball  (R=5). If you increase the radius, the relative error of this difference decreases. But IntegerPoints2  procedure is ill-suited for counting the number of points if the number of points is very large (the procedure seeking the points themselves).

This is a simple code, which finds the exact number of integer points inside the 6-dimensional ball of radius 20:

N:=0:

for i from -19 to 19 do

for j from -isqrt(400-i^2) to isqrt(400-i^2) do

for k from -isqrt(400-i^2-j^2) to isqrt(400-i^2-j^2) do

for l from -isqrt(400-i^2-j^2-k^2) to isqrt(400-i^2-j^2-k^2) do

for m from -isqrt(400-i^2-j^2-k^2-l^2) to isqrt(400-i^2-j^2-k^2-l^2) do

for n from -isqrt(400-i^2-j^2-k^2-l^2-m^2) to isqrt(400-i^2-j^2-k^2-l^2-m^2) do

if i^2+j^2+k^2+l^2+m^2+n^2<400 then N:=N+1 fi;

od: od: od: od: od: od:

N;  # The exact number of integer points in 6-D ball of radius 20

V:=(1/6)*Pi^3*20^6:  # The exact volume of 6-D ball of radius 20

evalf((V-N)/V*100);   # The relative error is 0.3 percent

                                                      329715877

                                                     0.3077223011

Another variant of Preben's example:

u:=4^(1/2)+sqrt(x^2)/(1+x);

map(`^`, u, 2);

u^~2;  # Maple just ignores  ~  operator

We see that  map  operator has a wide range of applications than  ~ operator.

@Carl Love   Thanks for the explanation.

@Carl Love   What is the meaning of this assignment  r:= table() ?

The next version works similar:

twinprimes:= proc(a::realcons, Lim::realcons)

local

     a0:= `if`(a::prime, a, nextprime(a)),

     b:= nextprime(a0), n, r;

          for n while b <= Lim do

          if b-a0 = 2 then r[n]:= [a0,b] end if;

          (a0,b):= (b,nextprime(a0))

     end do;

     convert(r, list)

end proc:

The procedure  IntegerPoints1  solves the same problem  as  IntegerPoints , has the same input parameters and the same output. But it is much faster and more reliable than  IntegerPoints .   IntegerPoints1  is based on brute force method  and doesn't use  LinearMultivariateSystem  or  isolve .

IntegerPoints1 := proc (SN::{list, set}, Var::list)

local SN1, sn, n, i, p, q, xl, xr, Xl, Xr, X, T, k, t, S;

uses combinat;

SN1 := SN;

for sn in SN1 do if type(sn, `<`) then SN1 := subs(sn = (`<=`(op(sn))), SN1) fi; od;

 n := nops(Var);

for i to n do

p := simplex[minimize](Var[i], SN1); q := simplex[maximize](Var[i], SN1);

if p = {} or q = {} then return {} else

if p = NULL or q = NULL then error "The region should be bounded" else

xl[i] := eval(Var[i], p);

xr[i] := eval(Var[i], q) fi; fi; od;

Xl := ceil~(convert(xl, list)); Xr := floor~(convert(xr, list));

X := [seq([$ Xl[i] .. Xr[i]], i = 1 .. n)];

T := cartprod(X); k := 0;

while not T[finished] do

t := T[nextvalue]();

if convert(eval(SN, Var=~t), `and`) then

k := k+1; S[k] := t fi; od;

S := convert(S, set);

if type(S, set(list)) then S else {} fi;

end proc:

Markiyan's example :

with(PolyhedralSets): with(ExampleSets):

rs := RandomSolid(6, 4):

#  the solid  rs  is increased 5 times 

IntegerPoints1(eval(Relations(rs), Coordinates(rs)=~(`/`~(Coordinates(rs), 5))), Coordinates(rs));

nops(%);

Thanks  Markiyan  for helpful criticisms and  vv  for the constructive suggestions that have been used in the procedure.

IntergerPoints1.mw

@vv  Thank you for your constructive suggestions. Using them  I'll try to rewrite the procedure  IntegerPoints  without   LinearMultivariateSystem  command.

@Markiyan Hirnyk  Thank you for your interest. The reason for the error is  the built-in command  SolveTools[Inequality][LinearMultivariateSystem] which was used in the procedure  IntegerPoints . See the simple example with the same bug:

SolveTools[Inequality][LinearMultivariateSystem]([x[1] > 0, x[2] > 0, x[1]+x[2] < 3], [x[1], x[2]]);

    Error, (in Utilities:-SimpleAnd) invalid input: a string/name list is expected for sort method `lexorder`

If you do not use indexed variables, all right:

SolveTools[Inequality][LinearMultivariateSystem]([x1 > 0, x2 > 0, x1+x2 < 3], [x1, x2]);

                             {[{0 < x1, x1 < 3}, {0 < x2, x2 < 3-x1}]}

@litun  If you are interested in only one real root  in the appropriate range  for each  A , use the  fsolve  command.

Details in  Roots1.mw

@vv  

restart;

rsolve(f(n+1)=f(n)+c, f(n));

                                                          f(0)-c+c*(n+1)

@mskalsi  Perhaps this is a bug. We see that the result depends on name of a variable:

simplify(diff(u(x,t), x,t), {diff(u(x,t), x) = 0});

simplify(diff(u(x,y), x,y), {diff(u(x,y), x) = 0});

                                  diff(u(x,t),t,x)

                                           0