Did you give a look to the IterativeMaps package?

What does your phrase "So we want to find a substitution that removes the time dependence from u. One way is to find the maximum and see how it moves. Here, the first solution gives what we want." mean precisely?

What kind of transformation are you talking about?

It is obvious that eq17 also writes

u(x, y, z, t) = 2*diff(L(x, y, z, t)), x$2);

# with 

L(x, y, z, t) = log(f, x, y, z, t)

So if you don't want u to depend on t you have to take f  to be for instance

f(x, y, z, t) = A(x, y, z)*B(t)

Is that what you need?

In such case have a look to  hfz_mmcdara.mw

If not forget it.

For_your_information.mw shows it worked well with Maple 2015: likely a code regression

@Kitonum 

"The variant ... more universal"

I agree but it wasn't the question.

BTW an answer coser to what the OP asked would have been

A := n -> Matrix(n, (i, j) -> a_||i||j);
 #or
A := n -> Matrix(n, (i, j) -> cat(a_, i, j))

Have a good day

@Kitonum 

Why didn't you write simply

A:=n->Matrix(n, symbol=a):

instead of 

A:=n->Matrix(n,(i,j)->a[i,j]):

?

Indeed the two definitions below seem equivalent

A:=n->Matrix(n,(i,j)->a[i,j]):
lprint(A(2))
Matrix(2, 2, {(1, 1) = a[1, 1], (1, 2) = a[1, 2], (2, 1) = a[2, 1], (2, 2) = a[2, 2]}, datatype = anything, storage = rectangular, order = Fortran_order, shape = [])

B:=n->Matrix(n, symbol=a):
lprint(B(2))
Matrix(2, 2, {(1, 1) = a[1, 1], (1, 2) = a[1, 2], (2, 1) = a[2, 1], (2, 2) = a[2, 2]}, datatype = anything, storage = rectangular, order = Fortran_order, shape = [])

@salim-barzani 

Remove the "&" character in the name of your file for it can be downloaded

and solve(%, R) returns 0 with multiplicity 2.
Explanation pde_mmcdara.mw

To complete @gkokovidis answer: as 0 is neither negative nor positve, it's worth looking to what happens when you used the Heaviside function:

restart
f := x -> Heaviside(-x)
x -> Heaviside(-x)
f~([-1, 0, 1])
                       [1, undefined, 0]
NumericEventHandler(invalid_operation = `Heaviside/EventHandler`(value_at_zero = 1)):
f~([-1, 0, 1])
                           [1, 1, 0]
NumericEventHandler(invalid_operation = `Heaviside/EventHandler`(value_at_zero = 1/2)):
f~([-1, 0, 1])
                           [   1   ]
                           [1, -, 0]
                           [   2   ]
NumericEventHandler(invalid_operation = `Heaviside/EventHandler`(value_at_zero = 0)):
f~([-1, 0, 1])
                           [1, 0, 0]
_EnvUseHeavisideAsUnitStep := true: 
f~([-1, 0, 1])
                           [1, 1, 0]

# Watch out: running _EnvUseHeavisideAsUnitStep := true:  
# makes `Heaviside/EventHandler`(value_at_zero = 0) ineffective
NumericEventHandler(invalid_operation = `Heaviside/EventHandler`(value_at_zero = 0)):
f~([-1, 0, 1])
                           [1, 1, 0]

Here are two tricks I use to ease the coding and avoid forgetting some backslashes (I used the substitute "&" in the second trick but you use any other character).
Based on @acer's answer:Sugggestions.mw

@acer 

I was completely off the mark in thinking that the solid line came from the way the procedure was called.
Now that you mention it, it seems to me that I've already noticed this effect due to the thickness of the line.

Thanks again for your comment and explanation.

@acer Thanks a lot.

I have two questions:

1. As no Warning, `i` is implicitly declared local ...l is issued while i is not assigned to something (Maple 2015), why did you declare i as local? As a rule I explicitely declare as local all the names for which such a warning is fired: does that mean that the procedures I write this way are not safe?
    
2. Could you explain me what happens in the las case

   ecdf(S, 'color'="green", 'thickness'=5)

   Why are the extreme segments in plain linestyle and not in dashed linestyle? I see that is related to the way you write the optional parameters but I still don't understand what happens here.

@dharr 

I have done a few attempts with "my" solution (assuming it  answer the OP's need, which I'm not sure):
 

initial_conditions := phi(0) = phi__0, D(phi)(0) = varphi__0;
sys  := {ODE, initial_conditions}:

Digits := 10:
sol := dsolve(
         sys
         , numeric
         , method=rkf45
         , events=[[phi(x)=0, halt]]
         , parameters=[phi__0, varphi__0]
       );


sol(parameters=[ two_numeric_values]);

interface(warnlevel=0): sol(1000):

XMAX := 1000:
xmax := sol(eventfired=[1])[];
if xmax=0. then xmax := XMAX end if:


DocumentTools:-Tabulate(
  [
    plots:-odeplot(sol, [x, phi(x)], x=0..xmax)
    ,
    plots:-odeplot(sol, [phi(x), v(phi(x))], x=0..xmax)
  ]
  , width=60
)

I was never capable to get plots closed to the the picture the OP gives.


Another test is based on a visual approximation of phi(x) from curve A:

test := 0.022*(1-tanh((x-150)/40))/2: 
plot(test, x=0..250);             # looks like curve A "phi(x) vs x"
plot([test, v(test), x=0..250])   # far from curve A "v(phi(x)) vs phi(x)"

Finally : v(Phi) is a decreasing negative valued function of Phi from Phi varying from 0 to 0.022 (ranges from curve A):

v(Phi);
plot(v(Phi), Phi=0..0.022)

So I can't figure out where the" Curve A pattern v(phi(x)) vs phi(x)" could come from?
 

@mary120 

Could you provide a file containing the values of x and phi(x) for curve A ?

Thanks

@salim-barzani 

In your initial qiuestion I simply followed step-by-step the paper you attached and it was clear and self consistent.
It's no longer the case here and I don't see the connection between this new paper and your attached worksheets.

For instance (file a[12]-pde.mw): you have now a nonlinear pde with integral terms whose unknown depends on 3 independent variables (a linear pde in x and t in the initial question).

• How can you get rid of the nonlinear terms?
• Do you want me to apply what I did in my answer to the linear part alone of this new equation?
• What are the coefficients you want to equate to 0?
• ... and so on.
   

Here is what you must take as my final contribution to this thread xxx.mw.
The procedure F generalizes procedure f in my initial answer, procedure `w'/w` has been update to handle more general situations.
I hope you will get more help for someone else.

By the way: Do you realize that you don't ask questions about Maple but that you ask helpers to do what is done in a paper, which requires a huge lot of time and enough knowledge of the underlying Physics? For you all this Physics is likely evidence but it's not the case for other Mapleprimes members, including me at first place.

 

@dharr 

... here is another point:

• You surely noticed this "strange" IC

  initial_conditions := phi(0) = 0.022, D(phi)(0) = 0.0;

  and probably wondered why the OP wrote this curious D(phi)(0) = 0 extra condition?
  (Which is why you use in your worksheet 

  sol := dsolve([ODE2, initial_conditions[1]], phi(x), numeric, range = 0 .. 10, method = rkf45, output=listprocedure);

  )

• I guess the OP doesn't try to solve ODE but diff(ODE, x): look to the patterns this could give:
   

  >  restart; >  # Define constants mu := 0.01: nu := 1 - mu: beta := 0.05: alpha := 0.3: M := sqrt(0.704): gamma1 := 0.001: # Define the function v(phi) v := f -> (1-alpha)*M^2 - (1-alpha)*M*sqrt(M^2 - 2*f) + mu*(mu + beta*nu)*(1 - exp(f/(mu + beta*nu))) +             (nu/beta)*(mu + beta*nu)*(1 - exp(beta*f/(mu + beta*nu))) + (alpha/gamma1)*(1 - exp(-gamma1*f)): # Define the ODE (dphi/dx)^2 + v(phi) = 0 ODE := diff(phi(x), x)^2 + v(phi(x)) = 0: ODE := diff(ODE, x): # Initial conditions: dphi/dx = 0 at x = 0, phi = 0 at x = 0 initial_conditions := phi(0) = 0.022, D(phi)(0)=0.0007;   (1) >  sys  := {ODE, initial_conditions}: xmax := 100: sol  := dsolve(sys, numeric, range = 0..xmax, method = rkf45): plots:-odeplot(sol, [phi(x), v(phi(x))], x=0..xmax); >  plots:-odeplot(sol, [x, phi(x)], x=0..xmax); >

 

Download rk4_mmcdara.mw

 

Finally; are we sure the OP doesn't want to solve 

ODE := diff(phi(x), x$2) + v(phi(x));

?
In its initial problem v(phi(0)) is Imaginary, and then phi(x) is a complex valued function (whose numeric value along a given branch can be ontained using the dsolve option complex=true).
At no time the OP suggests the phi(x) is imaginary (it even writes "phi and x are real"), which makes me doubt about the correctness of the term diff(phi(x), x)^2