Use multiple assignments as in the example below:

a, b := solve(x^2-1=0);

a;

b;

Another way:

plots:-display( plot3d(exp(x*y), x= 0..1, y= 0..1),   plots:-intersectplot(y=exp(x), z=0, x= 0..1,

y=0..exp(1), z=0..exp(1), color= red, thickness= 2, axes=normal));

restart;

p := proc()

local A, B, D, m, n, k, r;

A := Matrix([[0, 1, 0, 0, 0], [1, 0, 1, 0, 0], [0, 1, 0, 1, 0], [0, 0, 1, 0, 1], [0, 0, 0, 1, 0]]);

  r := LinearAlgebra:-RowDimension(A);

  D := Matrix(r); 

   for k from 1 to r-1 do

     B := A^k; 

       for m from 1 to r do

           for n from 1 to r do

               if m <> n and B[m,n] <> 0 and D[m,n] = 0 then D[m, n] := k;

               end if;

           end do;

        end do;

   end do;

   D;

end proc:

p();

Tools->Options->Interface->Default zoom

Download P.mw

In Maple you can use  seq  command.

Example:

<seq(1..10)>;  # This is a column vector

 or

<($ 1..10)>;

May be  GraphTheory[GetVertexPositions]  and   GraphTheory[Edges]  commands help you. The first command returns the coordinates of all vertices of a graph, and the second one returns the set of edges.

Another way:

x1 := solve(ln(x+2) = -1):   x2 := solve(ln(x+2) = 2):

plots[implicitplot](x*(y-ln(x+2)) >= 0, x = x1 .. x2, y = -1 .. 2, coloring = [blue, white], filledregions, gridrefine = 5, scaling = constrained);

A more automated solution can be get if you use  IntegerPoints  or  IntegerPoints1  procedure from here

eq1 := J+10*A+50*T = 500:

eq2 := J+A+T = 100:

IntegerPoints1({eq1, eq2, A >= 0, J >= 0, T >= 0}, [J, A, T]);

                                               [[60, 39, 1]]

RungeKutta4 procedure solves your problem for equation y'=f(t,y)  with initial condition  y(t0)=y0 .  The procedure was written on the basis of the formulas from wiki . Formal parameters: f  is the function  (t,y)->f(t,y)  defined by the right side of the equation, t0  and  y0  are defined by the initial condition, Tc  is the step,  T is the value of  the variable t  that t[N-1]<T<=t[N] ,  t[N] is the final value of  t , for which the procedure computes  y(t) . The procrdure returnes the list  [[t[0],y[0]], [t[1],y[1]], ... ,  [t[N],y[N]]] .

restart;

RungeKutta4:=proc(f, t0, y0, Tc, T)

local  N, t, y, h, L, n, k1, k2, k3, k4;

N:=ceil((T-t0)/Tc);

t[0]:=evalf(t0);  y[0]:=evalf(y0);  h:=evalf(Tc);  L[1]:=[t[0], y[0]];

for n from 0 to N-1 do

k1:=f(t[n],y[n]);

k2:=f(t[n]+h/2,y[n]+h/2*k1);

k3:=f(t[n]+h/2,y[n]+h/2*k2);

k4:=f(t[n]+h,y[n]+h*k3);

y[n+1]:=y[n]+h/6*(k1+2*k2+2*k3+k4);

t[n+1]:=t[n]+h;

L[n+2]:=[t[n+1],y[n+1]];

od;

convert(L, list);

end proc:

An model example: to solve the equation  y'=sin(t), y(0)=1  (the exact solution  y=2-cos(x) )

RungeKutta4((t, y)->sin(t), 0, 1, 0.1, 5);

plot(%, color=red, thickness=2, scaling=constrained, view=[0..5,0..3], labels=[t,y(t)]);

 Edit. The code has been corrected according to acer's comment.

First, we make the change  g(n)=sqrt(f(n))+1 , and easy to obtain a new recurrence relation

g(n+1)=g(n)+1  with the initial condition g(1)=sqrt(c)+1  (we assume that  f(1)=с ):

restart;

rsolve({g(n+1)=g(n)+1, g(1)=sqrt(c)+1}, g(n));

solve(%=sqrt(f(n))+1, f(n));

factor(subs({sqrt(c)=d,c=d^2},%));

f(n)=subs(d=sqrt(c),%);  # The final result

Carl, unfortunately your approach works only with integer degrees. Here is another solution that works with any degree:

P := 1+sin(e+e^1.5+e^2.5+e^3):

applyrule(e^n::numeric = `if`(n >= 2, 0, e^n), P);

                                        1+sin(e+e^1.5)

plots[implicitplot]((y-4)*(y-x^2) <= 0, x = 0 .. 2, y = 0 .. 4, filledregions = true, scaling = constrained, gridrefine = 5);

Addition 1. For details see help on  plots[implicitplot]  command  with option  filledregions=true 

Addition 2. Another general method you can see  here  (Picture  procedure).

This is only possible if  B  is not yet calculated:

B := Int(A(x), x = -infinity .. infinity);

op(1,B);

Standard way is using of  allvalues  command:

allvalues(SOLL[2]);