@Math-dashti It seems DEplot isn't set up for this. I'm not sure about fieldplot. Of course you could convert back to cartesian coordinates and use DEplot. odeplot works but is tedious.

Download polarphase.mw

@Math-dashti The output from PDEtools:-ConservedCurrents({de1,de2}, [u(t),v(t)]); correctly says any function of it will do. I stripped off the arbitrary function and just multiplied by that constant just to agree with @Rouben Rostamian's result. If something is conserved, then any function of it is also conserved.

@Math-dashti FYI, Maple can calculate the conserved quantity from the odes by

PDEtools:-ConservedCurrents({de1,de2}, [u(t),v(t)]);
-(1/2)*op(1, rhs(op(%)));

(The -1/2 just to convert to @Rouben Rostamian's form.)

@Rouben Rostamian  I'm fascinated by symbolic solutions, which was my main interest. Otherwise, I'm not convinced it is that useful. There have been many requests about it on Mapleprimes, so it was more a coding challenge.

@Alfred_F Here's one way. You always have to play around with these things.

Download test.mw

@Rouben Rostamian  If the point is that this is not useful as a numerical method, then I agree. But it does allow one to deduce the correct symbolic solution at x=0, which can then be verified and is well behaved. That presumably is a hint for the solution for general x.

In the post I wasn't enthusiastic about it as a numerical method, and the plot there shows divergence at relatively small x and t values for that example, though not as bad as here. Some of the literature I've seen only gives such comparisons for very short times, before the solution becomes poor.

Download LAD_Barenblatt0.mw

@Rouben Rostamian  The discontinuity may pose problems. Apparently "a given initial condition is consistent only with a specific set of boundary conditions"[1]; not sure if that applies here. Anyway, the method gives the same answer as a series expansion of the exact solution in t, so the algorithm appears to work. The hard part is deducing a simple closed form from the series solution, and I didn't make much progress here (I tried not to cheat since I knew the answer.)

LAD_Barenblatt.mw

[1] A. Garcıa-Olivares, Kybernetes, 32 (2003) 548.

@janhardo 

Yes, it can definitely be extended to other types of PDEs. Nonlinear parts that are more complicated than polynomials, such as exp(U*U[x]), might work with the existing code just by removing the polynom check, but I haven't tested it so I'm being careful. Additional time-dependent derivative terms can probably be added to L, and mixed time/spatial terms can probably be dealt with; these just need more code to detect and deal with the initial conditions. ODEs can be dealt with just by treating the highest-order term similarly to the time-dependent term here.

One difficulty is finding reliable test cases in the literature; the literature in this area is full of typos and missing information, for example figures that give values for most but not all of the parameters.

@salim-barzani The value you set for infolevel[LAD] determines how much extra information is printed out durung the calculation. You do not have to set it at all, in which case there is no extra information printed out. Or you can set it to 1, 2, 3 or 4 to include the information described at the top of the worksheet. Setting to 0 gives no information and setting to > 4 is the same as setting to 4. Many Maple commands, e.g., solve, will print information if you set infolevel.
 

In the code I used the userinfo commend to actually print out the information.

@KIRAN SAJJAN You didn't provide a worksheet. Did you change Sb? But if there are remaining errors they probably arise from other errors in the paper, rather than the Maple.

@KIRAN SAJJAN When I look at the worksheet it looks like H = Sb*D(Phi)(0)/(Le*Pr) instead of H(0) = Sb*D(Phi)(0)/(Le*Pr), but it is more than that.

Actually I think the original boundary condition is correct but the sign of Sb is wong. If I set Sb = -1 instead of Sb = 1, then I get correct looking plots. Later in the paper there is a plot where H is either positive or negative at eta = 0 for different signs of Sb.

Download dsolve2.mw

@KIRAN SAJJAN In the pdf, eq 13 is

and in your worksheet you have H(0) = -Sb*D(Phi)(0)/(Le*Pr). Now from the figures D(Phi)(0) is negative, and Sb, Le and Pr are positive, so H(0) will be positive, as it is in your worksheet. But according to the figure, H(0) must be negative. So the paper and your worksheet are wrong. I am suggesting you try removing the negative sign in this boundary condition and see if that works. Of course there may be other errors because of other errors in the paper or your worksheet related to the physics, but as far as I can see Maple is correctly doing what you asked it to do.

@KIRAN SAJJAN Your sahin_base_paper.mw won't run for me - in one of the boundary conditions (the one mentioned in my answer) you have "phis" instead of phi. The plots in that worksheet still have positive H at eta = 0 and you did not change the sign in the boundary condition. lt appears that you didn't read the comment at the end of my worksheet.

@mmcdara Thanks. I used to use approxsoln fairly frequently when I was solving these sorts of Navier-Stokes related systems, and it was always hit-and-miss. I think I was lucky here since the actual solution for H at eta=0 was only a small positive value not too different (relatively) from the negative value in the approxsoln. One can provide a digitized version if as in this case you know the solution, but usually that is not the case. Yes, ability to specify just one of the functions would be a nice feature.

@janhardo I agree that this is part of it. Large output should be suppressed because it takes time to decide it is too large to display.