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

@salim-barzani You omitted the eval(  , lambda=0) part, though when I added that it was still too slow.

@salim-barzani Here's how to do the plot. I don't know how to export the table; I usually generate data in Maple and write it out to a file for use in some other software. Hope someone else can help you.

 table-e2.mw

@janhardo Thanks for the exact solution and source. I'll work some more on comparisons.

@salim-barzani You have to recognize the general term manually. For simple cases guessgf can help (see here for an example).

Download Ex1.mw

@janhardo My procedure should work for cases with conjugates, which is what I think you mean by complex. My only problem is that I don't have any test cases with conjugates that have known exact solutions so I can check it.