@Preben Alsholm 

I could not reproduce some your steps. Maybe, the version of Maple or platform matter?

1. When I save the solution to *m file (Maple internal formft), read operator reads nothing from the file. Saving into a file with any other extension don't produce any issues (all works).
2. Syntax subs(S1,z(x)) is ideed simpler, I shall use it for further.
3. There is no problem with printing by showstat(`dsolve/numeric/hermite`) (as you metioned above). But assigning `dsolve/numeric/hermite`:=proc (................... doesn't solve the problem. After entering F1(0.12345)  the following error occurs: Error, (in F1) invalid input: diff received .13, which is not valid for its 2nd argument.

The cause is that the upper limit for variable y is not constant.

@acer 

Since 2D Output is a character style, I suppose, that alignment is defined by the following:

@Kitonum 

You are right, but for some reasons I would like to force Maple to work properly........

@Kitonum 

Thank you for advise, but, to my mind, the simplest way is the use of abacus istead of computer (it's a joke, don't be offended).

@ 

Thank you, now I'm aware of this approach.

@Carl Love 

Thank you, Carl. Indeed, Digits = 15 is a threshold for described issue. I would like to know, if it is possible to increase the precision of dsolve's solution in order to calculate integrals with higher precision (Digits > 15)?

@Carl Love 

Having tryed all ways above, I see, that the simplest method for calculation of integrals of powers of procedures (u^a) is to write own procedure of numerical integration based on some integration scheme with control of precision and convergence.

@Carl Love 

Here the result of second workaround:

It is a true curse!

@Carl Love 

18.02

@Carl Love 

Dear Carl, integration of u² does not work:

@Axel Vogt 

Thank you, Axel. Es funktioniert.

@Rouben Rostamian  

Thank you. Certainly, It’s a simple straightforward way and if it is not possible to calculate derivative of proc explicitly, it can be used.

@Preben Alsholm 

Of cause, dsolve returns the derivative of solution , I know this:

But I need to calculate the derivative of second order of solution ().

Is it possible?

@Carl Love 

Thank you for reply, I guessed that Matrices/Vectors are handled in special way in Maple sets.

Regards