Done with Maple 2015 (I verified that Maple 18 contains all the features for file 1 to work with this version).

If I interpret correctly what you say...

Cumulative_Histogram_with_Polygon_Line__Way_1.mw


In this Cumulative_Histogram_with_Polygon_Line_Way_2.mw worksheet the histogram is buit bars-by-bars using plottools:-rectangle instead of Statistics:-Histogram.
This avoids constructing the WeightedX list, which might be useful when Dados is a large sample, and offers more flexibility to produce smarter histigrams.

Illustration (note that the tickmarks could hve been displayed this way using the previous wotksheet).

So make your choice.


NOTE: the deprecated stats package (Maple 2015 and beyond) contained a useful feature enabling to draw histograms of weighted: Cumulative_Histogram_with_Polygon_Line__Way_1_stats_pckg.mw


 

As you use these two terms in the same question I understand you mean "envelope in the signalprocessing sense".

There exists a simple analogic device named Envelope Detector which, I quote Wikipedia,  "takes a (relatively) high-frequency signal as input and outputs the envelope of the original signal".

It is quite easy to draw the differential equation which models the input-output relation for this device (in the attached file R_d represents the resistance named R in the Wiki article, C_d the capacitance named C in this article and r_d is an auxiliary resistance left to the diode and which Wiki does not account for).

This auxiliary edo is added to your eqn1 and eqn2 differential equations.
To manage the two "charge" and discharge" stages of the envelope detector device I use a discrete variable GT(t) whose value +1 corresponds to the charge of the cacitance and value -1 to its discharge.
The values (-1 and +1) are controlled by events in the dsolve function.

Because V₂(t) is rapidly oscillating it is necessary to refine the odeplot outputs, which imposes in return to use the dsolve option range. 

File demodulation_toy_problem.mw shows how the envelope detector works. 
The solution depends on the three parameters R_d , C_d and r_d which have to be adjusted. Each choice gives an envelope but some seem more "reasonable" than others, or at least seem to fit more than others a mind image of what an envelope should be.

File demodulation.mw applies this envelope detector to your problem.

You will likely observe a quite large display time (which goes with a large amount of memory used), but this is the price to pay to a fully resolved solution.
Which explains why I did not go up to t=6*T on my memory limited machine.

__________________________________________________________________________________


In those files the 3-parameter model  is a physical one.
A simpler 2-parameter physical model results in setting r_d to some arbitrary value, for instance r_d = 1 (in true envelope detector one has R_d >> r_d ).
One may also think to a 2-parameter model (no longer that physical) whose parameters would be the two time-constants (one for the forward diode and the other for the reverse diode).

With this latter model the time-constant of the forward diode must be small enough for the envelope to follow faithfully the rise of V₂(t).
If, as it seems it is your case, the successive maxima of V₂(t) have increasing values (for t/T >~  0.9), then the time-constant for the reverse diode should be chosen large enough for  the envelope does not present significant decreasing between two consecutive peaks.

__________________________________________________________________________________


Finally,  demodulation.mw solves a 3-ode system but an alternative which could be interesting is to keep solving your initial 2-ode system first, and next dsolve the Envelope(t) ode alone while telling dsolve that V₂(t) is now a known function (see known option and how to use it in dsolve/numeric help page).
This could mimic the post-processing stuff you mentioned.

Nevertheless not that concise as MMA

Download Maple_vs_MMA.mw

Try for instance

limit__FG1 := limit(FG1, epsilon = 0);

    exp(eta[1]) * exp(eta[2]) + exp(eta[1]) + exp(eta[2]) + 1

and have a look to Limit ( help(Limit) )

@C_R 

I was writting a quite detailed comment to your question but an unexpected box failure made me lost all of it.
So for short, because I'm going to bed right now, I send you a code with several rationl fraction fits of your data.
Jut tell me if you are interested in it and I will develop all of this tomorrow.

RationalFractionFit.mw

... and the answer remains the same: How can expect that solve will find the the solution of a non linear system of 17 equations in 8 unknowns?
Unless the rank of this system is 8 (and there is no guarantee solve will even succeed) his is a dead end.
So you will find in the attached file an alternative method which enables finding a solution "in the least squares sense" (something I already told you about a couple of times several months ago).

test-F-p_sand15.mw

As several members already told you the result is unevaluated in some "ancient" Maple versions.

Nevertheless evaluation can be obtained if you observe that the exponent 

2*sqrt(tau^2-1)/(sqrt(tau-1)*sqrt(tau+1))

is equal to 2 if tau > -1 and -2 otherwise and if you use proper simplification assumptions. 
Evaluation.mw

PS: I suppose someone smarter than me could achieve the same result in a simpler way.

try

Sol := unapply(sol, x):
eval(convert(IC, diff), Sol(x0));

Or in a row

eval(convert(IC, diff), unapply(sol, x)(x0))

Download A_way.mw

By the way, why are you so prudish that you refuse to call  this "other software" MMA ?

The  modulus of solnum is the simplest case to deal with.
Given it's theoritical range of variation [0, +oo) the effective range the grid size "captures" can go from 0 to 10⁴ with a [300, 300] grid.

So it is impossible to get a smart representation of the level curves unless you lessen this range of variation. One simple way to do this is to apply the log function (or some surd function) to the expression to plot.
Using a "compressor"  such as log, for instance, introduces extra issues because the quantity to plot must be strictly positive (whichis why I concentrated on |solnum| ).
Note that a "compressor" is a trick which enables to get a smart plot, but the level curves or the colored region do not give a faithful representation.

The main idea to have a good spacing between level curves is to aggregate pieces of |solnum| contourplorts instead of |solnum| ina row.
Examples are provided in the attached file.

To reduce the size of the grids and still get a nice picture I suggest you to account for the z-periodicity of solnum (see the attached file).

At the end you will be able to get something like this by tuning the regularizer (f), the ranges of |solnum| pieces (r1, r2, r3), the contour lists or the coloring.


The plots are removed in the the attached worksheet (size 66 Mb instead)
plot-help_sand15.mw

This is not because Maple is powerfull that not using what  one does by hand is a useless strategy. Which in your case means that replacing k by a complex expression before trying to simplify this latter is a very bad way to proceed.

BTW: A syntax like  f(x) := something is not accepted in Maple 2015 (the version I use), So I used these f := x -> something notations instead. Feel free to modify them in your own code.

More of this there is no need to define Phi(x, t) the way you did because Phi(x, t) = C*Omega(x, t) where C is some constant

Download pdetest2_sand15.mw

After the M, b := ... line insert these to verify what b is:

b;


# So b[3] is not assigned:

b[3];
Error, Vector index out of range

The error comes from you writting

M, b := GenerateMatrix(eqs, [lambda[1], lambda[2]])

Which redefines b as a zero vector of dimension 2: thus  b = <0 | 0> and the b[3] term in M generates an error.

Cure: write instead

M, B := GenerateMatrix(eqs, [lambda[1], lambda[2]])

to keep b unchanged

S_sand15.mw

@nm 

As far as I can remember, exportplot, called from cmaple, has never worked properly for versions 15 through 2021, whether on Windows or Mac OS.

I always used this kind of workaround:

d := currentdir()

# f := cat(d, "/test.png »);  # doesn't work with maple 2015 (driver unknown)

f := cat(d, "/test.jpeg");

p:=plot(x^3, x = -8 .. 8, color = "blue",axis=[gridlines=[10,color="red"]]):

g := seq(plot(y, x=-8..8, color=red), y=[seq](-500..500, 100)), seq(plot([[x, -500], [x, 500]], color=red), x=[seq](-8..8, 2)):

pg := plots:-display(p, g):    

plottools:-exportplot(f, pg);  

Here is the file exported from Cmaple (Max OSX)

Note that the rendering still remains quite bad.

... which is probably non optimal and whose efficiency could likely be enhanced by someone who have enough time to spend on it.
Note that an obvious way to fasten this code is to add the numpoints option in the plot command. For instance numpoints=30 makes the code about seven times faster.

The code is generic in the sense it will enable you to

• plot any of your 6 functions
• versus any of one of the 4 parameters,
• giving to any one of the 3 remaining ones a list of values,
• and to the last 2 parameters a single value.

Generic_2d_plot_sand15.mw

The plots below are those you asked for

But as I said earlier the code enables you to draw any expression built from traceable quantities (those  listprocedure contains):


Several examples are given in the attached file, among them:

@Teep 

From Optimization/General/Options help page:

depthlimit = posint
Set the maximum depth of the branch-and-bound tree. This option is only used by the Optimization[LPSolve] command for integer programs. 
The default is the maximum of 2 and ceil(3*n/2).
The internal workspace allocated by the integer programming solver depends on the value of depthlimit.  If, during the computation, the workspace is found to be insufficient, an error is issued indicating that depthlimit must be increased.

The attached worksheet presents some tracks on a simplified problem (5 jobs and 5 machines):

1. For this n=26 variable problem the default value of depthlimit is 114 and wasn't large enougth to avoid the

   Error, (in Optimization:-LPSolve) maximal depth, 500, of branch and bound search is too small; use depthlimit option to increase depth

2. Using depthlimit=1000 returns a solution with performances

   memory used=3.47MiB, alloc change=1.93MiB, cpu time=3.70m, real time=3.71m, gc time=0ns

3. Using depthlimit=600 returns the same solution with better performances

   memory used=2.21MiB, alloc change=1.62MiB, cpu time=5.10m, real time=5.13m, gc time=0ns
   Which suggests to use only the necessary depthlimit value, no more (which is rethorical)

4. Using depthlimit=500 returns the same error as with the default depthlimit value

5. Using again the same depthlimit=500 value but restricting the number of nodes to explore (nodelimit=10^5)
   gives almost the same solution than depthlimit=600 with nodelimit=0 ... but it's  about seven times faster:
   memory used=1.97MiB, alloc change=1.55MiB, cpu time=43.61s, real time=43.67s, gc time=0ns

So I suggest you not to keep increase depthlimit (which penalizes the computational time), but set it to  a value reasonably larger than the default one (182) and increase progressively nodelimit until the  optimal solution no longer evolves.

5 jobs and 5 machines  Some_tracks_5_5.mw   !!! read red text below
5 jobs and 8 machines  Some_tracks_5_8.mw   !!! read red text below

When opening again those same files with Maple 2015 (the version I use) suprious errors appear due to some invisible character in some command blocks.
Here is the (partially) rewritten code for the 5 jobs and 8 machines case 
Some_tracks_5_8_(new).mw

XT_sand15.mw provides a very simple solution.

Nevertheless changing the x range or the tvalue will very likely make this error to appear

Error, (in dsolve/numeric/bvp) Newton iteration is not converging

Despite it is a very classical error for which some recipes do exist to fix it, see dsolve/numeric_bvp/advanced help page and the so-called continuation method, it can be hard to find the correct way to use this trick.
I blindly tried some possibilities but even I can't say any of them gave satisfaction.
Another possibility might be ti initoiate the dsolve process by giving some approximate solution  (search for approsoln in dsolve/numeric/BVP help page).

Generaly this kind of issue cannot be fixed by applying blindly those recpes and the one which is best placed to use them in an enlightened manner is the one who understands the physics underlying the differential system to be solved, i.e., you.