I am long harassed by finding a numerical method to the problem of rotating disk.

Here, the name "rotating disk" problem refer to a set of four ordinary differential equations modelled from the axis-symmetric flow from a rotating disk (see P163 ~ 168 of the book <<Viscous Fluid Flow>> by Frank M. White):

The easiest way to solve this equation is to write the equations in finite difference form and Newton iterate. Otherwise try the book "Two point boundary value problems" by H.B. Keller which I vaguely recall discusses the related Falkner-Skan problem. You could also try Schlicting or any other book on boundary-layers (your problem is usually referred to as the von Karman similarity solution).

Yes, I have consulted that book. However, in Falkner-Skan problem, there is only one value to be shooted, so that's relatively easy.

I even solved Blasius equation numerically only by linear shooting.

But for this rotating disk problem, two values are to be determined.

Schlicting's classical book also didn't talk of how to numerically get F'(0) and G'(0), it mentioned that those two valued were obtained by the method of series expansion and match. But today this method is clearly two complicated and outdated. With the power of modern computer, it should not be that difficult to solve it numerically.

Actally, I have tried to mail Prof Keller about this problem, coz he is an established expert in boundary-value problem of ODE; but i got no reply.

Have you tried my ealier suggestion of Newton iteration on the finite difference version of the equations - there is no need to know either F'(0) or G'(0) in this case?

I think you're being a bit unfair on the series solution method - have a look in either the Journal of Fluid Mechanics or the European Journal of Mechanics B/Fluids to see that these methods/ideas are still used. Also have a look at