according to my knowledge., by using Group theoretic method analysis, or lie group methods.  The system of governing nonlinear partial differential equations with boundary conditions reduces to a system of ordinary differential equations with appropriate boundary conditions.

Note that your equation for h_2 is of first order in r, hence much easier to solve. It seems that assuming R time independent is consistent with the equations, and all conditions. Hence, to solve the equation for h_1, I would make the ansatz h_1(t,r) = t^(-1/4) f_1(x), where x= (r/R)^2. That should lead you to a solution f_1(x) = C_1*(1-x)^(1/3), where C is a constant (probably the one you found in the textbook, whose solution would be wrong if R depended on time -- I think). 

Next make the ansatz h_2(t,r) = t^(-1/4) f_2(x). That gives you a first order equation in x for f_2(x).

And I honestly don't think you will get any closer to the solution if you start by studying 3 semesters of group theoretic methods ;-)

Note added: If you assume that h_2(t,r) = K h_1(t,r), for some constant K, then your equations are essentially identical. You only have to determine the constant K.

Hi Kim, I guess there is a printing error in the second PDE. The parabolic term is surely the second derivative of h2, right?

If yes, my suggestion: forget about analytical solutions. Invest your time better into numerical solution and do it not yourself, seek professionals as partners. Try first the stationary case and introduce a source instead of the r.h.s. zeros. My guess, you have a flow problem in a tube or so? Besides, are you allowed to lift the secret of the physical nature of h1 and h2?

I guess, by using symmetry method, you can reduce the order of your PDE by one. Then you can solve by choosing appropriate ansatz. Or you can use an equivalence transformation to make your PDE integrable. This will be handy so far I understand. 

my suggestion is lie symmetry method for solving this equation and you can easily applied to lie symmetry method with the aid of maple package program.

@Kim: I repeat, your second PDE is surely written incorrect. The second term there should contain a second derivative of h2. But this printing error aside, you have to give us the initial conds if you really wish to get help. Another issue. If your PDEs are conservation laws then please tell us which sort of symmetry you have, cylinder/axial or spherical symmetry? In both cases I cannot see that your system is complete. And your diffusivity is increasing from r=0 upwards? This is a turbulent problem or why is diffusivity increasing? More info  needed. Strange problem.

@Kim: I make a U turn as  I now realized that you have full cylinder symmetry and as an initial condition something like a delta/dirac function? Besides the wrong written parabolic term in the second of your PDEs the initial condition needs to be specified. The allegedly analytical solution (k) exhibits a singularity at t -> 0. What is still irritating is the solution-dependent diffusivities of h1 and h2. They are not identical so that turbulent mixing may be excluded. Yes, this is a hard task and will take time to solve it. Besides, I assume that h2 exhibits a similar structure in terms of separation of the variables. There shouldn't be too much difference. So from a practical point of view I fully agree with Kare Olaussen. Don't invest too much into Lie groups etc.. But please verify that in the second PDE the term dh1/dr should be  replaced by dh2/dr.

Unless I have made a mistake, expression (k) is NOT a solution of PDE (j).  The parabolic shape may be some form of leading order approximate solution, in which case I derive the time dependence should be t^{-1/3}. So you need to check where these PDEs come from and re-derive the results.

It would appear that the first PDE (j) is a thin fluid film approximation, with circular symmetry,  strong normal gravity, and negligible surface tension.  The second PDE (l) looks like some sort of linearisation of (j) as if one were seeking solutions h(r,t,eps)=h_1(r,t)+eps*h_2(r,t)+...  but if so then I believe there is missing term from PDE (l).  You need to check its derivation also.

Apart from the question of whether equation (l) is correct or not, there is the more serious problem that equations (j) and (l) are not dimensionally correct (unless the notation is extremely misleading). It is really obvious when you look at the equations, but I did not discover it before evaluating the normalization constant for h_1, which cannot be made to agree with the solution (k).

I'm really sorry for my late reply. Thank you all for replies.

@Helmut Ziegfeld Baumert

you mean there are 2 boundary conditions so second derivative of h2?

these pictures are my systems. that 2 immiscible liquid layers spread radially

and in my guess(though I don't know how textbook authors solve it), mass balance can be used as an initial condition

@A. J. Roberts

I take your point that the time dependence should be t^{-1/3}. but I guess using mass balance can convert order of t

@Kåre Olaussen

no dimension problems in (j) and (l), v(nu) means kinematic viscosity[m^2/s] h(height)g(gravity acceleration constant) in these eqautions.

@Hany A. Hosham ,  Abdus Sattar Mia, Sait San

Oh there's some theories to learn.

@Kim: Thanks for this clarification. Tonight I am going northward for three weeks on an island where there is (hopefully!!!) NO internet so that I cannot longer follow your interesting problem. Good luck & crossing fingers!

Mj> v(nu) means kinematic viscosity

OK, I interpreted v as velocity.

I already told you how to find h_2 from h_1. But since there has been raised some doubts that h_1 is correct (it is) I upload a complete analysis.

Note added: I must say that an incredible amount of irrelevant and misleading advice was posted to your question.

@    Helmut Ziegfeld Baumert

Thanks for the help you have given!

@Kare

Brian G Higgins gives information about h1

Although I haven't fully flipped through it. I want to share this with you

This notes are made by Brian G Higgins

Mj> information about h1

Looks correct.

Regarding the similarity solution: It is likely that some of the assumptions use to derive its equations break down for very small  t, where this solution becomes singular. Probably also for very large t.

Hi.in this gard refer to P.Brandimarte book

Regarding the comment "It is likely that some of the assumptions use to derive its equations break down for very small t, where this solution becomes singular." True.  But for a range of similarity problems one can prove that the similarity solution emerges from quite general initial conditions---including ones that are within the domain of validity of the derivation of the equations.  That is, often the similarity solution is widely relevant irrespective of its notional singularity at time zero. 

For one example of this see the paper   S. A. Suslov and A. J. Roberts. Similarity, attraction and initial conditions in an example of nonlinear diffusion. J. Austral. Math. Soc. B, 40(E):E1–E26, Oct. 1998. http://journal.austms.org.au/ojs/index.php/ANZIAMJ/article/view/453.