one possibility

restart;
F:=x*y^2-x;
G:=x*sin(Pi*x);
solve([F,G],[x,y]);

_EnvAllSolutions:=true;
solve([F,G],[x,y])

or

restart;
F:=x*y^2-x;
G:=x*sin(Pi*x);
solve([F,G],[x,y],allsolutions=true);

I would not do this. This is the wrong way to program.  

restart;
sol:=solve(x^2+3*x+4=0,x);
assign(('A','B')=sol):
         

A better way is to just assign all the solutions to one variable. A list. Then use indexing to access each solution.

The less variables there are, the more control on your code you have., You do not want a function with 100 variables in it.

The above assumes there are 2 solutions ofcourse. In code, if you want to do this in a more robust way, you'd need to first find the number of solutions, then come up with, at run-time, variable names on the fly, and do the above, so they match.

For example:

Assign solution to list and access sol using [] indexing or iteration

restart;
sol:=[solve(x^2+3*x+4=0,x)];
for N,item in sol do
    print("solution ",N," is ",item);
od;

#or
sol[1];  

sol[2];  #etc..

or

restart;
sol:=[solve(x^2+3*x+4=0,x)];
for i from 1 to numelems(sol) do
    print(sol[i]);
od;

Only one variable is needed. sol to store all solutions. No need to come up with other variable names.

I found the problem. I did not look too close before.

To get the arrows, you need to have coupling between the two ODE's. Before you had 

    OnDE:=diff(On(t),t)=(rhs(solEr))/500*60-On(t)/1500*60;

where you replaced the actual Er(t) above. This removed the coupling between the two ODE's.

Instead, you need to just write Er(t) in the second ODE, and Maple will figure it out.  Now Maple see the two ODE's make a system.

If you do not have coupling, then these were 2 separate ODE's and that is why the arrows were missing.

restart;
ErDE:=diff(Er(t),t)=30-Er(t)/500*60;
OnDE:=diff(On(t),t)=Er(t)/500*60-On(t)/1500*60; 
initvals := [[Er(0)=0,On(0)=0], [Er(0)=500,On(0)=0],[Er(0)=500,On(0)=1500]];
DEtools:-DEplot( {ErDE,OnDE} , {Er(t),On(t)}, t=0..150, initvals, numpoints=100, linecolour=[blue,red,yellow]); 

if I understand right, see if this works for you. To make it go by 1,3,5,7,... just add "2" after, so it jump by 2 each time, starting from 1.

restart;
h:=x->sin(x)/x;
A:=n->Int(mul(h(x/a),a=1..2*n+1,2),x=0..infinity);

A(2);
A(3);
A(4);

You can add value() to evaluate the integral at end.

in numerical integration, you can't have parameter with no numerical value in the integrand. (how is Maple going to find numerical value if a has no value?)

Assign a some value before calling int, then it will work

a:=.4;
f := (x, y) ->(1/2)*a*(sinh(y-x^2)+tanh(x-y^3)); 
sol := evalf(int(int(f(x, y), x = -6 .. 5), y = -5 .. 5));

or

restart;

Explore( 
   evalf(int(int((1/2)*a*(sinh(y-x^2)+tanh(x-y^3)), x = -6 .. 5), y = -5 .. 5)),
   a=1..10,initialvalues=[a=1],size=[800,60]
);

You can also use Explore, and change the raduis using the slider with the mouse.

restart;
Explore(plots:-display( plottools:-sphere([0,0,0],r),view=[-3 .. 3,-3..3,-3..3]),
        r=0.1..3.0,initialvalues=[r=0.5]):

always use exact numbers, unless there is a reason not to

a := 12/100;
kn := n -> (n + 1/2)*Pi/a;
g := (n, x) -> cos(kn(n)*x);
g(1, a);

   0

It makes it easier if you replace  by say A before.

restart;	
s1 := Ls*cos(omega*t + phi__l + theta)*omega + sin(omega*t + phi__l + theta)*Rs;
s2 := sqrt((omega*Ls)^2+Rs^2)*sin(omega*t + phi__l + theta + arctan(omega*Ls/Rs));
s1:=subs(omega*t + phi__l + theta=A,s1):
s2:=algsubs(omega*t + phi__l + theta=A,s2):	

And now

s2:=simplify(expand(s2)) assuming positive:
evalb(s2=s1)

                true

You can also replace the above by

s2:=simplify(expand(s2),symbolic);
evalb(s2=s1)

            true

another option, if you want to do it more manualy

restart;
roll_col_1:=rand(10..20):
roll_col_2:=rand(50..100):
col_1:= Vector([seq(roll_col_1(),i=1..6)]);
col_2:= Vector([seq(roll_col_2(),i=1..6)]);
A:=<col_1|col_2>

You need to be careful which convention to use for spherical coordinates. There is the Physics one and there is the mathematics one. See https://en.wikipedia.org/wiki/Spherical_coordinate_system

Using the Physics one (more common)

evalf(int(int(int(ln(r^2 + theta^2 + phi^2 +1),r=0..2),theta=-Pi/2..Pi/2),phi=0..2*Pi))

97.2256142300

Using the math convention for the angles

evalf(int(int(int(ln(r^2 + theta^2 + phi^2 +1),r=0..2),theta=0..2*Pi),phi=-Pi/2..Pi/2))

97.2256142300

if you use 0..Pi instead of -Pi/2..Pi/2 you get

evalf(int(int(int(ln(r^2 + theta^2 + phi^2 +1),r=0..2),theta=0..2*Pi),phi=0..Pi))

105.5237634000

Not sure which convention you want. Wikepida page above describes these more.

ps. only numerical integration worked in Maple. 

restart;
ode:=diff(y(x),x)+k*piecewise(y(x)>=0,y(x),y(x)<0,0) =sin(x);
dsolve(ode)

Can't solve it analytically. You can solve this nuemrically only

restart;
k:=1;
ode:=diff(y(x),x)+k*piecewise(y(x)>=0,y(x),y(x)<0,0) =sin(x);
sol:=dsolve([ode,y(0)=1],numeric)

Maple used to have many of these cases, but it got better over years handling them.

A trick to bypass this issue, is to delay the substitution of the actual function until after you obtained the series solution. As follows.

Using a generic f(x,y) first, then this function is replaced by the actual value of the function after the series solution is obtained. This makes life much easier for Maple pdsolve.
 

Download pdsolve.mw

Looking at Wikipeida, it looks like you used different formula. The formula with the origin at center is

You wrote

Where the eccentricity you used is 0.25.

Fixing this, and increasing e, you get ellipse.

e:=0.25; b:=1;
p1:=plot(b/sqrt(1- (e*cos(theta))^2), theta = 0 .. 2*Pi, coords = polar, scaling = constrained,color=red):
e:=0.75;
p2:=plot(b/sqrt(1- (e*cos(theta))^2), theta = 0 .. 2*Pi, coords = polar, scaling = constrained,color=blue):

plots:-display([p1,p2])

Try this. Need to use := and not = for assignment. And why write 1.  instead of just 1 ? This is not Matlab. If you want real numbers, you could always convert the final answer to float, and it is more accurate this way.

restart;

H := proc(n::integer)::real; 
   local s, i; 
   s := 0; 
   for i by 2 to n do 
       if i <= n - 1 then 
          s := s + 1/(i*(i + 1)); 
       else 
          s := s + 1/i; 
       end if; 
   end do; 
   return s;
end proc;

A := proc(n::integer, T::procedure)::real; 
       local d1, d2, s; 
       d1 := T(n + 1) - T(n); 
       d2 := T(n + 2) - T(n); 
       s := T(n + 1) - d1*d2/(d2 - d1); 
       return s; 
end proc;

I have not used spacecurve before in Maple, but looking at help, I think you need to do this.

restart;
f:=(x,y)->2*x^2+3*y^2-x*y-6;
subs([x=1,y=1],diff(f(x,y),x));
subs([x=1,y=1],diff(f(x,y),y));
p1:=plot3d(f(x,y),style=patchcontour):
p2:=plots:-spacecurve([x,1,f(x,1)],x=-10..10):
p3:=plots:-spacecurve([1,y,f(1,y)],y=-10..10):
plots:-display([p1,p2,p3])