This is not the error but a warning:

restart:

with(plots):

Digits:=35:

ini1:=D(x)(0)=0, x(0)=1:

dsys:=diff(x(t),t,t)+(x(t)-2)*diff(x(t),t)+5*x(t)=0;

dsol1 :=dsolve({dsys,ini1}, numeric, abserr=1e-9, relerr=1e-8, maxfun=0);

plots:-odeplot(dsol1, [x(t), diff(x(t),t)], t=0..2, axes=boxed, color=black, linestyle=1, tickmarks=[6, 6], axes=boxed, titlefont=[SYMBOL,12]);   # The plot  x'(t) against  x(t)

The issue is the appearance of multiple roots. See

solve((x+y+z)^2=0, z);

                                   -x-y,  -x-y

A workaround:

plot3d(solve((x+y+z)^2=0, z, dropmultiplicity = true), x=-100..100, y=-100..100);

A:=(cos(sqrt(2)*t)-cos(omega*t))/(omega^2-2);

``(1/denom(A))*numer(A);

For better visibility animation done in 2D. It is easy to rewrite in 3D by using circular symmetry about Ox-axis.

A := plot([-1, 1], x = -10 .. 10, filled, color = "Yellow"):

B := plots[animate](plot, [[2*r*cos(phi)/Pi, 2*arcsin(sin(r*sin(phi)))/Pi, phi = 0 .. 2*Pi], thickness = 3, color = blue, scaling = constrained], r = 0 .. 5*Pi, frames = 100):

plots[display](A, B);

 Addition - animation in 3D:

A := plot3d([x,cos(t),sin(t)], x = -10 .. 10, t=0..2*Pi, style=surface, color = "Yellow", transparency=0.75):

B := plots[animate](plot3d, [[2*r*cos(phi)/Pi, 2*arcsin(sin(r*sin(phi)))/Pi*cos(t), 2*arcsin(sin(r*sin(phi)))/Pi*sin(t)],phi = 0 .. 2*Pi, t=0..2*Pi, style=surface,color = "Niagara Azure", scaling=constrained, numpoints=10000], r = 0 .. 5*Pi, frames = 100):

plots[display](A, B, axes=none);

Better quality of plotting without additional option  grid  enables the use of spherical coordinates.

Compare:

restart;

A:=plots[implicitplot3d](x^2+y^2+z^2-2^2, x = -2 .. 2, y = -2 .. 2, z = -2 .. 2, style = surface, axes = normal, view = [(-2.7 .. 2.7)$3],  grid = [30, 30, 30]):

B:=plot3d([2*cos(s)*sin(t), 2*sin(s)*sin(t), 2*cos(t)], s=0..2*Pi, t=0..Pi, style=surface):

plots[display](<A | B>, axes = normal, view = [(-2.7 .. 2.7)$3]);

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.