@Mariusz Iwaniuk I got the following message during the calculation:

"FindRoot::lstol: The line search decreased the step size to within tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the merit function. You may need more than 15.` digits of working precision to meet these tolerances."

Haven't you got such a message?

@Markiyan Hirnyk thank you, I have not tried the Roots-function before. My recent experiences with this function for higher order zeroes (n>60) is that it is slower than Next->Zero. Moreover, sometimes it gives the same zero many times. Probably, this is because the value of f(x,k) is too close to zero for large values of n. Certainly, f(x,k) can be rescaled appropriately for larger values of x to circumvent this issue.  

@Mariusz Iwaniuk perfect, thanks! 

By the way, here is the data of Fig.1 of the original post:

Beyond these data I got the instabilities. I hope that I can reproduced this matrix with both of your algorithms (in Maple and Mathematica) and see which one is faster and more reliable.

I will also try to extend it to a table of 40x40 or more, although my progamming skills are not sufficiently good such that I could understand your algorithms in any detail. 

@Markiyan Hirnyk typo correction: replace k^2 by k. (sorry)

k := 5; t := time(); L := [6.543504327];

for j to 10 do L := [op(L), RootFinding:-NextZero(proc (x) options operator, arrow;

exp(x)*HeunC(2*x, 1/2, -1/2, -x^2, k*(k+1)+1/8, .999999999) end proc, L[j])] end do; L;

[6.543504327, 8.100773700, 9.742070947, 11.44667626, 13.19955834, 

  14.98981792, 16.80947993, 18.65263039, 20.51481776, 

  22.39264021, 24.28345989]
 

This is corresponding to line 6 in the matrix of my original post.

@Mariusz Iwaniuk I have applied the Maple function "NextZero" with

                                    f(x,k) := exp(x)*HeunC(2*x, 1/2, -1/2, -x^2, k * (k+1) + 1/8, .999999999)

for k=0,1,2,.... 

The reference is https://arxiv.org/abs/1807.09128 equation (25).

EDIT: k^2 replaced by k in the above function.