@JoyDivisionMan

Here is an alternative way to what I named the "pointwise construction" of the self-convolution of function f.
Based on the SignalProcessing:-Convolution function this computation should be much faster than the  "pointwise construction" I presented in my answer (generally numeric convolution and correlation use the FFT algorithm), but requires some carefulness to correctly rescale the result.

Here is the code (X is the same 'X' than in my reply with step=0.01).

a := Array(f~(X), 'datatype' = 'float'[ 8 ] ):
b := a:
c := SignalProcessing:-Convolution(a, b):

N := numelems(c):
K := 2*(max-min)(X)/(N-1):

#----------------------- raw convolution
# dataplot(c, style=line);

plots:-display(
  #----------------------- rescaling to get 'c' in true values
  plottools:-transform(
    (x, y) -> [x*K, y*K])
    (dataplot(c, style=line, color=red, legend="SignalProcessing:-Convolution")
  )
  ,
  plot([seq([x, Exact(x)], x in X)], legend="exact", color=blue)
  , view=[0..4, default]
)

Result

@JoyDivisionMan 

Detailed explanation using a toy problem

Download Detailed_explanation_with_a_toy_example.mw

"This system of nonlinear equations was proposed by an artificial intelligence named Alice." means absolutely nothing out of context.

Alice's answers are nonsensical out of context and the last "Practical applications" item is one of the moset stupid thing I ever read.
By the way Alice forgot to mention this set of equations could also have a practical application in predicting 

She should have stayed  in Wonderland where she was more brilliant.

@Andiguys 

Question_1_sand15.mw and reply_1_sand15.mw give indeed different results because they do not solve the same equations.
To be clearer the equations you solve in your initial question and the comment to my answer are not the same.

Pay attention to what you did and read First_and_second_problem_compared.mw carefully.

@Andiguys 

Read carefully the attached file

Download reply_1_sand15.mw

@nm 

Thanks for sharing your experience

Maybe a track: the mw file I wanted to upload contains debugger outputs: if I remove them the file can be uploaded.

For information: I get no error using Maple 2015, even for a little bit complex display.

With_Maple2015.mw

@acer 

You're right that A^N/2 * B^N/2 by (A*B)^N/2  are not always equivalent and I had even thought to write "Assuming that A^N/2 * B^N/2 by (A*B)^N/2  holds..."
But I was also guided by the expressions of the two terms the OP wanted to prove the equality, and without this assumption I was pretty sure the proof could not be established
So I thought that this equivalence goes without saying.

But again you're right, I sould have been detailed in my answer.

@FDS 

Is there a policy to mark your latest post as best?
Yes, you should see a trophy at the top right of my answer: just click on it.

See you next time!

@FDS 

The updated worksheet  Stress_Strain_Curve_2.mw
And an illustraton of what you can do with it in the spirit of what you show in your pptx

About polynomial regression... there is really no mistery behind it.
The first question is: Are you familiar with the "multiple" linear regression in the case of Q regressors X₁, ...X_Q?
If it is so, "regressing" a dense polynomial of degree Q-1 with a single indeterminate Z, is nothing but using the "multiple" linear regression while setting X_n = Z^(n-1), n=1..Q+1.
To give an example, suppose you want to fit a dense quadratic polynomial a+bZ+cZ²: simply introduce auxiliary variables X₁ = 1, X₂ =Z, X₃ =Z² and regress your dependent quantity no longer on  a+bZ+cZ²:  but on aX₁+bX₂+cX₃...and that's it!

If you are not familiar with "multiple" regression I advice you to read any elementary textbook/course/mooc... on the subject.

Beyond that an important question often appears: what is the degree of the polynomial to use in order that the fitted polynomial fulfills two contradictory requirements:

1. Fidelity to the data: otherwise said the fitted polynomial must be "close" to the data or, mathematically speaking, the residual sum of squares between the data and their restitution by the fitted polynomial must be small.
   Obviously, as a dense polynomial of degree D-1 depends on D parameters, a set od data of size D sill be exactly fitted by this dense polynomial.
   We use to say here that the model is over-parameterized
2. Stability (or robustness, or generalization capability...): let us take this previous sample of size D and this same polynomial of degree D-1.
   This polynomial may have D-1 real roots and if, unhappily, some of them are located within the range of the regressor (Z), then the fitted polynomial will present undesired overshoots in between consecutive values of the regressor.
   In fact it is extremely likely that the fitted polynomial will behave this way.
   If you want to avoid this phenomenon use a polynomial of degree 0 or 1.

So, to satisfy the fidelity to the data requirement you must take a large degree, close to the sample size, but if you are more concern by "stability" or "smooth fit", then you must take low degree polynomials.
There exist a lot of litterature about strategies which balance these two contradictory requirements. I already gave you a few Wiki references on the subject.

What I wrote in my worksheet is that each branch is a sample of size about one thousand.
So a dense polynomial of degree 5, for instance, is very far ffrom being an  over-parameterized polynomial. This balance stuff I spoke abour above matters only when the number of parameter of the model (not necessarily a polynomial one) is of the order of the sample size, let us say 20% or 25% of this latter.
You would be right considering model selection (here degree selection) if you were investigating polynomials of degree 200 for example. But it is very farfrom being the case.

As I wrote in this last worksheet a mor important criterion is that two polynomial models intersect within L5 and EP5 ranges if you do not want to create unealistic approximated cycles.

@FDS 

Feel free to tell me if this goes to the right direction
Stress_Strain_Curve_2.mw

I understand the notion of cyclic loading-unloading , but you want to assess the "surface area" of each cycle, and so you need to define more precisely what characterizes a cycle in terms of geometry: where does it start (likely the intersection point between a charge path and the next/previous discharge path) and where it ends.
So I propose you my interpretation for the first cycle only... I am not going to go further if I am not on the right track in me doing this.

My Maple version is too old to manage workbooks. Could you simply upload the stress-strain data?

@janhardo 

Beautiful, I vote up

Error messages are pretty clear: you have to specify an approximate solution.

For more detais execute help(dsolve,numeric_bvp,advanced) and help(dsolve,numeric,BVP) in a worksheet and search for the word approxsoln.

Defining an approximate solution can be quite tricky and the best way to do this is to construct a simplified version your problem using a simplified physics (do a dimention analysis, identify second order effects, discard them), to solve it, and to use this solution as approxsoln (which can be either a procedure or an array).
Your are likely the better placed to do this because it is you who knows the physics behind your equations.

For instance, what is order of magnitude of
compared to other terms? Can you simplify it?

Or can you approximate BCs like?
Can you simplify them from a physical point of view?
As a rule a BC on the first derivative which nonlinearly depends on the second derivative makes a BC problem complicated to solve.

My advice is: begin to do this by yourself and come back to us if you get Maple difficulties to sove your problem.

@JoyDivisionMan 

By the way, the attached worksheet explains in more details how to get the pdf of Y = cos(X) where X  is a uniform random variable with support (-𝜋, +𝜋).
Two approaches are used: a direct construction of PDF(Y) and an indirect one by derivation of CDF(Y).

To compute the cdf of, let's say, X = sin(X) (which is obviously the same as Y) without doing all the stuff contained in the attached file (solve/allsolutions, branch identification, aggregation of partial pdf/cdf over the different branches, ...) it is much simpler to observe that sin(X) = cos(𝜋/2 - X) then define an auxiliary random variable A = 𝜋/2 - X and replace everywhere it occurs  f_X(something) by f_A(some other thing) which is equal to f_X(𝜋/2 - something)

Download cosine.mw