@acer 

Thanks, I will remember that

@love-algebra 

Probably not the simplest way to do this, but it works...

restart:
expr := 3*L*r/m/l/s+L/m:
MyRule := beta=L/m;

subs(L=solve(MyRule, L), expr);

eval(L[1, 1](r, theta), u__0(r, theta) = f[1, 1](r, theta));

# or

subs(u__0(r, theta) = f[1, 1](r, theta), L[1, 1](r, theta))

@Thomas Richard 

Answered from home, 

Thank you Thomas.

" Since your optimization call succeeds both with method=ego and method=diffevol, there is probably no need to try other methods and options - apart from curiosity, which is okay. ;-)" 
Right, but the curiosity wasn't accidental: In my real application the solution of the ODE system depends on 10 parameters
(with quite different ranges of variation, but I normalized both of them between 0 and 1 to avoid numerical difficulties).
The first time I used GOT with diffevol I got an error message concerning the time limit.  Not the message you find in GetLastSolution help page, but an error message (I will be able to send it tomorrow) wich, from memory, said "error in convert/radical ... in < here some procedure, not always the same > ".
So I began to try other things.
Finally the problem was solved by setting the option "localrefinement" (sorry, I cite from memory) to false and by using diffevol ("ego" is much slower, so the idea to reduce the size of the design of experiments to build the response surfaces)

If you are interested I could post the true application tomorrow.

@zaidoun 

You're right, sol_z (output (5)) is the solution.

Please: note that the solution I obtained here was already found by Kitonum.

@zaidoun 

The command 

substitutes any occurences of z by exp(u), thus u becomes the "new variable".
The idea, following the initial remark Carl made, is to reduce the search range which extended to 5*10^10 to somethong smaller (z = 5*10^10 means u = log(5*10^10) ~25).
One problem is that your initial range for the imaginary part of z is -1000.. 5*10^10.
To cover this range with the change of variable z --> exp(u) is impossible. So, what I suggest you is:

1. solve F(z) for a reasonable range: previous answers showed RootFinding works well for re=-1000..1000 and im=-1000..5000.
2. Next solve G(u) on the extentendedrange  re=-1000..1000, im=log(5000)..log(5*10^10)

The command 

enables recovering z solutions sol_z from u solutions sol_u.
Remain sol_u is a list, so the command exp(sol_u) won,'t work. To recover sol_u you need:

1. to build a loop over all the elements of sol_u and compute the exponential of each of them,
2. or to apply exp to all the elements of sol_u in a single shot:
   1. either by using the tilda symbol (see the corresponding help page for the tilda symbol)
   2. or by writing sol_z := map(exp, sol_u)

Finally the two last command evaluate G(u) for each u in sol_u and F(z) for each z in sol_z.
And YES, the result is (then 4 times the same, a thing that ought to be verified)
 

@zaidoun 

Maybe you could find some ideas here


 

Download TryThis.mw

Don't have Maple 2018 right now, here is the result obtained with Maple 2015.2

 

Download pde-2015.mw

@Carl Love 

It's the kind of thing you can do when there's nothing more interesting around here

• I use to visit the site https://oeis.org/?language=english and I went to search if the sequence
  1, 28, 110, 384, 1316, 4496, 15352, 52416, 178960, 611008 ...  would not have interesting properties.
  Unfortunately this sequence is not identified, but, if you set f[1]=1 and f[2]=M, then f[n] = A[n]*M+B[n] where
  A = 0, 1, 4, 14, 48, 164, 560 ... and B = 1, 0, -2, -8, -28, -96, -328 ...
  Sequence A[n], from n> 1, is dubbed A007070
  Sequence B[n], from n> 2, is dubbed A106731
  They are both related and many characterizations for both of them, including generating functions and expansions of rational fractions are given.
  Unfortunately I wasn't capable to find any clue to prove the sequence generated with M=28 doesn't contain any square (except f[1] of course).
  Maybe you would be more successful?
   
• So I had a little fun, here are a few results I got
   

  Download fun.mw
  
  Observations: 

  • For f[1]=1 only 25 values of f[2] in the range [2..675] contains (at least) one f[n] which is a square when n varies from 1 to 1000

  • The corresponding value of n is rather small... maybe it is even unique?

  • One also generate sequences with f[1]=p and f[2]=m (m>p). One finds that for some couples (m, p) there exist different values of n such that f[n] is a square

@Carl Love 

Thanks Carl.
I agree, these alse prompts "are a total nuisance to remove"

@phiost 

No problem.

In case you would be nevertheless interested in coloring "3D level curves" and identifying them by a legend you could find a few ideas here.

contourplot2.mw

@Carl Love 

Just a remark about the "prompt", 

If you copy-paste your first conditional sequence into a worksheet (worksheet mode) then its five line are preceeded with a prompt.
But, as they are in the same single execution group, the result is the one expected.
The prompts obtained by this copy-paste "are the same" as those you get when applying menu>edit>join on multiple groups.
Then it seems there are two kind of prompts, saying the "active" ones which begin a group and the "passive" ones inside a group.
So maybe it's notcompletely true to say "in other words, a single command prompt".
 

@vv 

Smart, as Iwrote to silvergasshoper I had done that years ago.
Your solution is definitely better, I thumb up.

A point of detail...

In the case where the thermal conductivity is discontinuous, your scheme doesn't guarantee the (physical) continuity of the thermal flux (unless the mesh size is adpated in a correct way left and right of the discontinuity points).
The classical method to avoid such unphysical solutions consists in deriving finite difference approximations that explicitely account for the continuity equations for the temperature and the thermal flux.
One can show that the correct thermal continuity at its points of disconinuity is the harmonic mean of its values left and right to these points.

On coarse grids the "harmonic" method is by far the better. On finer grids its superiority on the classical method (the one you used) tends to lessen. More of this the "harmonic" method is more time consuming and may be found complex to implement efficiently on unstructured meshes.
 

Download HarmonicMean.mw

To complete Tomleslie's answer, you could be interested to know that the value of  time "t0" can also be found by using the events feature of dsolve
(type "events" in the help and look to the  dsolve, numeric, Events (events) help page)

With tomleslies' notations ("events" are supported by rkf45 [default}, ck45 and rosenbrock methods)

sol1:= dsolve( 
               { sys, x(0)=0, y(0)=0, z(0)=cos(1), w(0)=sin(1)},
               type=numeric, method=rosenbrock, events=[[[x(t)<0, And(-y(t)<-5.5, y(t) < 6.5)], halt]]
             );

LargeEnoughTime := 20;  # set LargeEnoughTime to a value larger than the one returned by the warning message
sol1(LargeEnoughTime): 

EndTime := op(sol1(eventfired=[1]));
sol1(EndTime);

EndTime := 7.84831340510192