@chrisc 

because I still can't reproduce, see the attached (where I have also inserted a couple of 'version-checking' commands)

oddOne.mw

Note that I am running Maple 2020.2 on 64-bit Windows7. Exactly which Maple version/subversion are you running, on which OS?

@chrisc 

attached. Use the big green uparrow in the Mapleprimes toolbar

@ik74 

see the attached worksheet

bigodes3.mw

which will produce the plots

@vv 

it states

find the integer numbers a, m, n so that the given equation has integer solutions?

I'm not convinced that the number of 'integer solutions' is specified: although I grant that the OP also specifies

"Can this equation has four integer solutions"

But, does this mean - are there four integer soltuion for a, m, n, x?

Or does it mean are there four integer solutions for 'x' when a , m, n are also integers?

I assumed one of these - and you assumed the other

@vv 

You may want to consider (for example)

{a = 21, m = 15616744, n = 73, x = 22}

for which the equation

a^2*x^2/(x-a)^2+x^2 = m/n:

returns

213928 = 213928

See the code below to see how to generate  an infinite number of such solutions

  restart:
  eq:=a^2*x^2/(x-a)^2+x^2 = m/n:
  sol:=simplify~(isolve(eq));
#
# Pick any two integers for _Z1, and _Z2,
# substitute these in 'sol' and then evaluate
# 'eq' for the resulting values of a, m, n, x
#
  eval(sol, [_Z1=22, _Z2=73]);
  eval(eq, %);

@ik74 

where the "extra" execution group plots xi[j](t), versus diff(xi[j](t),t). Again because there are only three distinct eqautions, ony three distinct curves will be produced

The attached file

bigodes2.mw

will produce the curves

@Muhammad Usman 

to your file Q2.mw which appears to do what you want

Q2ans.mw

To transfer a matrix from one worksheet to another, you will have to use the ExportMatrix() command in the sheet where the target matrix is generated, and then the ImportMatrix() command in the worksheet where the matrix is desired

@Muhammad Usman 

but maybe(?) one of the options in the attached

genMAt2.mw

@Scot Gould 

1. For simple problems, it doesn't make much difference, whether you use output=listprocedure, or (the default) output=procedurelist
2. Sometimes, depending on subsequent calculations, it can be more convenient to use output=listprocedure.
3. However (generally?) it is less "efficient" to use the output=listprocedure option. If one has a system with multiple dependent variables (or high order derivatives of a single dependent variable), output=procedurelist will essentially compute all available dependent variables simultaneously. Using output=procedurelist to compute (say) three dependent variables, means that these will be computed "sequentially" - and this takes longer

So for the following "toy" example

dsys:= { diff(x(t),t,t) = -2*lambda(t)*x(t),
         diff(y(t),t,t) = -2*lambda(t)*y(t)-Pi^2,
         x(t)^2+y(t)^2 = 1,
         x(0)=0, D(x)(0)=1/10, y(0)=-1, D(y)(0)=0
       }:
dsol1 := dsolve(dsys, numeric):
CodeTools:-Usage(plots:-odeplot(dsol1, [[t, x(t)],[t, y(t)],[t, lambda(t)]], t=0..10)):
dsol2 := dsolve(dsys, numeric, output=listprocedure):
xf:=eval(x(t), dsol2):
yf:=eval(y(t), dsol2):
lf:=eval(lambda(t), dsol2 ):
CodeTools:-Usage( plot( [xf(t), yf(t), lf(t)], t=0..10)):

the ttimings are

memory used=276.80KiB, alloc change=0 bytes, cpu time=16.00ms, real time=14.00ms, gc time=0ns
memory used=13.55MiB, alloc change=0 bytes, cpu time=218.00ms, real time=209.00ms, gc time=0ns

Note the cputime/realtime/memory increase for the second option, using output=listprocedure

For this example, it doesn't really matter very much because the overall times are so short. However for more complicated systems where "efficient" generation of such plots might take a few seconds, such a 15X speed hit would be very noticeable!

@666 jvbasha 

Some examples

1. Yout latest post specifies a condition on the dependent variable 'theta'. This variable occurs nowhere else in your problem. So presumably this condition is completely irrelevant
2. The final plot you show seems to be parameterised by the variable 'delta'. The variable 'delta' occurs nowhere else in your problem
3. The function h(z) is defined for sigma<=z<=sigma+1. What value does it take outside this range?
4. In order to produce the final plot, you must have numerical values for every variable in your problem, eg 'xi', 'eta[1]' (and possibly others). What values are you using?

If this problem comes from a textbook/paper whatever, it might be easier if you just scanned and posted the textbook/paper here. Use the big green up-arrow in the Mapleprimes toolbar

@JAMET

it outputs

[.5000000000, 1.118033988*I]
Warning, data could not be converted to float Matrix

immediately before the plot. The complex number obtained for Ip is the source of the warning message

See the attached

yesItdoes.mw

@minhthien2016 

to vv's original code - see highlighted line in the following

  restart;
  ti:=time[real]():
  XY:=seq(seq(`if`(irem(x*y,x+y)=0,[x,y],NULL), y=2..x-1),x=2..30):
  N:=0:
  for xy in XY do
      x,y:=op(xy):
      for a from 2 to 10 do
          for b from 2 to a-1 do
              r:=a/x  + b/y; m:=numer(r); n:=denom(r);
              if n<=10 and m>=2 and m<=n and igcd(m,n)=1 and igcd(a,b,m)=1 
              then N:=N+1;
                   sol[N]:=[x, y, x*y/(x + y), a, b, m, n]
              fi
          od:
      od:
   od:
   time=time[real]()-ti;
  'N'=N,entries(sol);

This reduces the numberr of solutions to 64 - as obtained by Kitonum (but it runs even faster, because the reduction in the loop index measn there are many fewer cases to check!

@DarkMath 

you state correctly that for,
# Start snip

r*e_r + 0*e_th + 0*e_ph

Now, the function converts this to cartesian coordinate system with [0,0,0] being the root. 

The conversion of this vector to cartesian coordinate system must be

[r*sin(th)*cos(ph), r*sin(th)*cos(ph), r*cos(th)]

# End snip

Now if I set th=0 in the above I get

[r*0*cos(ph), r*0*cos(ph), r*1]; #or
[ 0,          0,           r  ]

So what is the problem?

It is not too difficult to produce the sphere which is is tangent to the edges, rather than the faces.

The following code produces both such spheres

  restart:
  with(geom3d):
  _EnvXName := x:
  _EnvYName := y:
  _EnvZName := z:
#
# Define the "centre" of the tetrahedron
# Note that [0,0,0] is just a kind of obvious
# choice - feel free to pick any other value
#
  point(o,0,0,0):
#
# Define the "in-radius" of the tetrahedron
#
  inRad:=2:
#
# Define the tetrahedron with  the centre and
# in-radius specified above
#
  tetrahedron( tet, o, in_radius=inRad):
#
# Define the sphere with the centre and radius
# specified above
#
  sphere(sph0, [o, inRad]):
#
# Compute the sphere which is sangent to the edges
# of the tetrahedron
#
  sphere( sph1, [o, sqrt(add((add(vertices(tet)[1..2][])/2-~coordinates(o))^~2))]):
#
# Draw the tetrahedron, the in-sphere, and the sphere
# which si tangent to the tetrahedron's edges
#
  draw( [ tet(  color=yellow, transparency=0.5),
          sph0( color=green,  transparency=0.5),
          sph1( color=red, transparency=0.5)
        ]
      );
#
# Check the details of the tetrahedron and the
# two spheres
#
  detail(tet);
  detail(sph0);
  detail(sph1);

Generating (amongst other things!!) the figure below

And the issue of an a sphere which is tangent to every edge of a tetrahedorn is more difficult than a sphere which is tangent to every face of a tetrahedron, so feel free to completely ignore the above answer.

I will make a resolution to read questions more carefully