The Task model seems inappropriate for the OP's purpose, which is performing a computation three different, independent ways. I see the Task model as being fundamentally recursive. I've implemented Roman's idea in my Answer below, but in the basic Threads model rather than the Task model.

@Axel Vogt The title was my fault, and I've changed it now. I'd assumed that the OP's final line of output was what was being claimed as erroneous, as is usual with a bug report.

@Mariusz Iwaniuk I'm sorry. Your inappropriate use of evalf threw me off. I thought that you meant that the bug was that evalf didn't return 0 like the others. Usually when reporting a bug you should make the final line of output the one that shows the bug, unless some other line of output is explicitly pointed out.

Do you have any initial and/or boundary conditions?

The command evalf doesn't supercede value. After you apply evalf to an inert indefinite integral, you still have an inert indefinite integral. Note that the integral sign is grey, not blue. There is no bug.

@maple2015 You wrote:

Please download the attached file

I did. Your file is difficult to work with because everything is jammed into one execution group, so the output only shows at the end. Please separate it into execution groups with only one significant command per execution group. And please display the results after those significant commands by ending them with a semicolon rather than a colon.

I recommend that you set Digits to 15 until you get the bugs worked out. That way it'll run much faster. And also forget about timing it until the bugs are worked out.

The value of C[2] is assumed to be zero for relative error checking ( as you recommended)

I didn't exactly recommend that you neglect setting the absolute error tolerance. I said that if you were only going to use one of relative and absolute error tolerance, then use relative. When the scale of the problem is on the order of 1, then it makes sense to set both, and to the same value.

1- It is reasonable that by assigning a little value to C[1], the little tolerance of C[1] for little amount of absolute maximum step size in interval [0,1], should be ineffective to the answer. By taking the value C[1]=Float(1,-15) and C[7]=1E-3, the answer is  0.229491 but for C[1]=Float(1.-14) the answer is 0.221455. Is it show the instability of solution?

I don't fully understand lsode. I think that it's not really intended for dynamic use of the solution function. Rather, I think it's primarily intended to integrate the ODEs to a specific stop value. I recommend that you read this recent thread concerning lsode, paying particular attention to the Answer and Replies by Allan Wittkopf.

So, why do you want to use lsode rather than any other of the numerous numeric methods?

2- I use the Implode and Explode commands in line 14. Is there a better way for obtaining the expression from 'RootOf' without using the String Tools?

Yes, it's trivial with the command op:

EX:= op(%);

It's pretty rare that StringTools are the best way to pick apart algebraic output.

3- Can I solve the ode in line 12 directly by maple commands?

Yes, it's easy. It's not really an ODE because it contains only the first-order derivative and doesn't contain the primary function undifferentiated. Thus, you can use ordinary solve to solve for diff(u(y), y).

Indeed, the whole process of solving, interpolating, integrating, and plotting can be reduced to

ode:= U__y^n = -G*(1-y)/mu+tau[Y]*(1-exp(-a*U__y))/mu:
F__EX:= r-> fsolve(eval(ode, y= r), U__y= 0..infinity):
CodeTools:-Usage(
     plot(
          'Int(F__EX, 0..x, epsilon= 1e-5)', x= 0..1, numpoints= 11,
          color= violet, thickness= 2, title= "Laminar Flow Velocity",
          gridlines= false
     )
);

Here it is in a worksheet. Please compare this plot with the plot that you obtained by your means.

Download Laminar.mw

@spreka The purpose of fnormal is to convert numbers that are close to 0 to actually be 0. This makes the remaining numbers, the true eigenvalues, much easier to see in the Matrix.

@acer Thank you for the improvements to my hasty code. They certainly expand its scope of applicability.

@Thomas Richard I emphasize that those FAQs apply only if the error occurs when starting Maple. The vast majority of kernel connection errors are caused by kernel bugs, not by firewall or configuration issues. Indeed, I've never seen a case of the latter, nor seen a personal report of such a case in any public forum.

@Joe Riel Pausing the output with a simple readstat could be useful because it lets you read the output that has been generated so far. This could be especially useful if the program has several traces and infolevels set---my primary debugging tools---generating reams of output.

There are a great many possible reasons for that. You'll need to post your code. Usually it is caused by an improperly trapped error in the internal Maple code (the kernel code) that causes the death of the kernel.

You asked:

When we use the relative and absolute errors (C[1] and C[2]) simultaneously it leads to the more accuracy or using absolute error is sufficient?

There's no easy answer to that. There is a correspondence between the relative error and the number of significant digits. Specifically, a relative error of .5*10^(-d) corresponds to d significant digits. If the magnitude of a result is less than one, then the relative error tolerance is more constraining than the absolute error tolerance; if the magnitude is greater than one, then that is reversed. I'd recommend never setting the relative error tolerance greater than .5e-4. The setting of the absolute error tolerance must depend on the scale (magnitude) of the expected results. If you were going to use only one of the two, I'd recommend using the relative. (I welcome any arguments to the contrary.) Note, however, that when the result is 0, achieving a certain relative error can be difficult or impossible.

Why after solution C[12] to C[21]  (the components of array corresponding to output data) are zero?

I see what you mean, but I don't have any answer. I think that it's a bug. I think that lsode is not frequently used, so not much maintenance is done on it. If it were being maintained, then Arrays (instead of arrays) would've been allowed many, many years ago.

@spreka After the iteration, the eigenvalues are not "randomly placed around the diagonal." Their placement is highly structured, as this shows:

Download QR.mw

@spreka Again, please provide a small example of such a matrix C.

@acer (I'm sure that you know this already; this is just for other readers.) The referenced code from Stack Exchange has a severe bug. There, the iteration is equivalent to

while not stopping criterion do    
     (Q,R):= LinearAlgebra:-QRDecomposition(C);
     C:= Q^+.R.Q
end do;

However, in the OP code, the iteration is correct.