@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.

Please provide a small example of what you consider to be a "bidiagonal" matrix. I've heard of tridiagonal, but I don't think that bidiagonal is well-defined. If a matrix has only one diagonal other than the main diagonal, then its eigenvalues are obvious: They are just the main diagonal.

@Kitonum I believe that the OP is refering to something like this:

restart:
M1:=146996733613391:
M2:=1348471408813:
teks:= CodeTools:-Usage(numtheory[cfrac](M1/M2,  'quotients'));

memory used=23.62MiB, alloc change=22.90MiB, cpu time=31.00ms, real time=28.00ms

31 ms seems to me to be an outrageously long time for such a trivial computation.

@marc sancandi wrote:

PS : For Carl Love :  at x=L1 the OP has not written boundary conditions (as you said) but just Temperature and heat flux continuity conditions

Although it may not be a "boundary condition" in the sense of a true physical boundary or in your textbook senses, Maple is counting any condition that uses a specific value of x as a boundary condition. Three different values of x are used. That is the one and only reason for the error message. The lack of a Dirichlet boundary condition is irrelevant to Maple at this point in its analysis of the problem. That's all that I was saying.

Maple's numeric PDE solver is quite limited. You've specified boundary conditions at three x values: x=0, x=L1, x=L1+L2. The solver can only handle boundary conditions at at most two x values.

@one man I guess that you're refering to Draghilev's method. It would be great to have a general procedure for that in Maple.