Adaptive is spelled thus. That might help with the Googling. A search on "adaptive Runge Kutta" (without quotes) turns up several relevant articles that explain the theory and show some code in other languages. You may need to translate that to Maple.

@Markiyan Hirnyk Any three points are concyclic. So if the six centers degenerate to three, no further proof of concyclicy is required.

@Markiyan Hirnyk 

It is a degerate case because P1=P6, P2=P3, and P4=P5. There are only three centers, so they are concyclic.

It is a challenge to do the integral accurately for even a single value of r (other than r=0). Changing the value of digits, epsilon, or Digits makes big changes in the answer, indicating that the answers are garbage.

The integral in the other Question is easy because of the strong exponential decay.

I think that some transformation will be needed for this integrand. Also I think that the interval should be split, e.g., 0..1 and 1..infinity, and then the infinite part translated to a finite interval.

I am skeptical about the Mathematica plot because of the sharp point near r=2.

@sarra

Change the line

while (x[n])<b) &(h>hmin)  

to

while x[n] < b and h > min do

No parentheses are required, `and` must be spelled out, and the word do is required.

Change display to printf, and put its argument in double quotes rather than the single quotes that you have now.

@acer Does setting epsilon or digits have any effect when using a NAG method? Isn't everything done in hardware floating point? Yes, I know that theoretically it could reduce the number of evaluation points, but does it actually?

@sarra Now I need to know what you intend the return value of the procedure to be. As it is, the return value will be the last computed value of y.

Assigning to the parameter y is problematic. What did you intend by that?

@Markiyan Hirnyk The number option is not needed in this case, because the program can easily count that nops(DE13) = 2. If you do supply the option, then its value must be 2. I don't understand why that particular error message is issued, but setting number= 0 is ridiculous because it means that there are zero equations in the ODE system.

Clearly the OP had in mind something other than counting the number of equations by specifying number= 1000.

Quoting the help

number specifies the number of equations in the procedure-form first order system....

The number option is used to specify the number of equations, and is not needed when the equations are supplied algebraically (as you have them here). Obviously you had something else in mind by using number= 1000. Perhaps you meant the numpoints or numsteps option.

@sarra error is a reseved word in Maple. Replace data[error] with data[err].

Your values for b1, b2, and b3 are correct. b1 can be any positive number less than sqrt(72). b3 can be any number greater than sqrt(72).

If b=0, there is no damping, and the period is 2*Pi/sqrt(4.5). If b < sqrt(72), then the quasi-period is 8*Pi/sqrt(72-b^2). Note that this reduces to the former answer when b=0.

@teh_allchemist Your ig, the gearbox setting, is a discrete variable: its values are limited to {1,2,3,4,5}. Solve for each of the five values separately. I did it for ig = 5, and the solution only took 81 seconds. I guess the others would be about the same. And the solution was "only" 16 screens wide by 1 screen tall.

@Entvex Ah, I see. You want to extend the function periodically. That's a bit more complicated.

f:= t-> piecewise(0 <= t and t < 2, 2-t, 2 <= t and t <= 3, 2*t-4):

plot(f(t-3*floor(t/3)), t= 0..6);

@Mary Your only real mistake is that Array declarations need to be made with parentheses, not square brackets:

X:= Array(1..n);

You should also make i a local variable.

After that, your procedure works.

The %a format code is a powerful catch-all.