Download A_way_mmcdara.mw

Examples

Download with_plottools.mw

You can take inspiration from this to produce more general graphics, or build your own contour drawing procedure.

Mathematically a steady state of some ode diff(s(t), t) = F(t, s(t)) means that diff(s(t), t) = 0 far from some value t=T, in which case the steady state is F(T, S(T))=0.

A weaker version iconsists in sying thatthe solution s(t) of some ODE has reached a quasi steady state at time T, if |s(T+dt) - s(T)| < epsilon for some "small" positive epsilon 
If s(t) keeps oscillating far from some value t=T, you cannot say that s(t) is in a quasi steadystate... unless you use a value of epsilon od the order ot the oscillation amplitude... which seems quite unreasonable.

As you are in this latter case it's not reasonable (IMO) to say you have reached some steady state.

But maybe you mean "periodic regime" instea of "steady state regime"?
If it is so here is a simple way to visualize the transition to periodicity. The idea is that the spectrum of a periodic regime is discrete.  So if you use a moving window of a signal, compute its spectrum, and observe a finite nimber of rays (or almost) from some window shift T, than you may reasonably consider you have enter a periodic regime.

This idea is illustrated below. Nevertheless it needs to be refined by using a proper windowing (shape, length) and tuning some parameters too. 
So please take this as a first step only

2_mmcdara.mw

Given the results in this file one might say that a "from the point of view of the spectrum amplitude [and given the parameters I choosed] a steady state regime is reached from t ~ 350).

Note: the width of the two rays (in fact only one if one plots only one half of the dft) comes from the t-step (=1) I used. Better results are likekely to be expected iif one replaces 

tau=1+10*(first-1)..10*(first-1)+N

by

tau=1+10*(first-1)..10*(first-1)+N, step

where step < 1 (0.1,  0.01... feel free to try this)

Download anim.mw

(A) Build different curves from point z to point z0 and, NEXT, plot a closed curve arround (which can be quite simple)

(B) OR build the closed curve and AFTER draw several strictly internal curves (which is more complex).

Other questions:

1. Is that for illustrative purpose alone or is the closed curve is defined by some equation?
2. Do you want to do math or just drawing... this latter can be easily done using PowerPoint for instance ?
   Did you look to Insert/Graph in Maple's menubar?

To adress item (A) above you can use this procedure:

Download SeveralPath.mw

@MathGuru 

(apart an absolute value in equation (8) one could likely explain the origin)

A few checkings and a numeric evaluation with a random value of t shows that H is almost certainly the exponential matrix of I*t*B

Download MatrixExponential.mw

L := [9,15,21,25,27,33,35,39,45,49,51,55,57,63,65,69,75,77,81,85,87];

# C is a list of 2^numelems(L) elements, so the "code" below is not intended
# to long lists L.

C := combinat:-choose(L):

N := 97:
A := map(c -> if add(c) = N then c end if, C)
[[9, 15, 21, 25, 27], [25, 33, 39], [25, 27, 45], [21, 27, 49], 

  [15, 33, 49], [9, 39, 49], [21, 25, 51], [15, 27, 55], 

  [9, 33, 55], [15, 25, 57], [9, 25, 63]]

Longer lists require a more subtle code Ihave no idea of.

The package ScientificErrorAnalysis is perfectly adpated to do "First Order Sensitivity Analysis" (FOSA).
In FOSA one uses only a central value of the target expression (which you name R) and its first derivatives (just as you do).
But, contrary to your point of view where the parameters have only infinitesimal variations you never give values to, FOSA defines each parameter by a central (nominal, or reference, if you prefer) value and an indication about its range of variation.
From a statistic point of view one might say that each parameter is defined by its mean and its standard-deviation.

If your gradient approach and FOSA share a lot of elements, FOSA goes far beyond as it accounts for the ranges of variation of each parameter.
More of that ScientificErrorAnalysis enables accounting for possible correlations among parameters.

Here is an illustration of what you can/could do with this package

With_ScientificErrorAnalysis.mw

Factored expression find-parameter_2_mmcdara.mw

with Maple 2015 provided you don't load usless packages: Q3_mmcdara.mw

I guess you mean "mathematical functions" ?
In such case you could use
?initialfunctions

To get all these names:

C := FunctionAdvisor(function_classes):
AllMathFunctions := op~(FunctionAdvisor~(C)):  # generates several lprint

AllMathFunctions;
[cos, cot, csc, sec, sin, tan, cosh, coth, csch, sech, sinh, 

  tanh, arccos, arccot, arccsc, arcsec, arcsin, arctan, arccosh, 

  arccoth, arccsch, arcsech, arcsinh, arctanh, arccos, arccosh, 

  arccot, arccoth, arccsc, arccsch, arcsec, arcsech, arcsin, 

  arcsinh, arctan, arctanh, cos, cosh, cot, coth, csc, csch, exp, 

  ln, sec, sech, sin, sinh, tan, tanh, Beta, GAMMA, binomial, 

  doublefactorial, factorial, lnGAMMA, multinomial, pochhammer, 

  Psi, harmonic, lnGAMMA, KelvinBei, KelvinBer, KelvinHei, 

  KelvinHer, KelvinKei, KelvinKer, AiryAi, AiryBi, HankelH1, 

  HankelH2, AiryAi, AiryBi, BesselI, BesselJ, BesselK, BesselY, 

  HankelH1, HankelH2, KelvinBei, KelvinBer, KelvinHei, KelvinHer, 

  KelvinKei, KelvinKer, AiryAi, AiryBi, BesselI, BesselJ, 

  BesselK, BesselY, HankelH1, HankelH2, KelvinBei, KelvinBer, 

  KelvinHei, KelvinHer, KelvinKei, KelvinKer, cos, cosh, cot, 

  coth, csc, csch, sec, sech, sin, sinh, tan, tanh, ChebyshevT, 

  ChebyshevU, GegenbauerC, HermiteH, JacobiP, LaguerreL, 

  LegendreP, Chi, Ci, Ei, Li, Shi, Si, Ssi, FresnelC, FresnelS, 

  Fresnelf, Fresnelg, dawson, erf, erfc, erfi, KummerM, KummerU, 

  WhittakerM, WhittakerW, CylinderD, CylinderU, CylinderV, 

  CoulombF, CylinderD, CylinderU, CylinderV, Ei, FresnelC, 

  FresnelS, Fresnelf, Fresnelg, GAMMA, HermiteH, KummerM, 

  KummerU, LaguerreL, WhittakerM, WhittakerW, dawson, erf, erfc, 

  erfi, EllipticCE, EllipticCK, EllipticCPi, EllipticE, 

  EllipticF, EllipticK, EllipticModulus, EllipticNome, 

  EllipticPi, InverseJacobiAM, InverseJacobiSN, LegendreP, 

  LegendreQ, ChebyshevT, ChebyshevU, ChebyshevT, ChebyshevU, 

  EllipticCE, EllipticCK, EllipticE, EllipticK, GaussAGM, 

  GegenbauerC, JacobiP, LegendreP, LegendreQ, LerchPhi, 

  SphericalY, arccos, arccosh, arccot, arccoth, arccsc, arccsch, 

  arcsec, arcsech, arcsin, arcsinh, arctan, arctanh, ln, 

  LommelS1, LommelS2, AngerJ, LommelS1, LommelS2, StruveH, 

  StruveL, WeberE, AiryAi, AiryBi, AngerJ, BesselI, BesselJ, 

  BesselK, BesselY, ChebyshevT, ChebyshevU, Chi, Ci, CoulombF, 

  CylinderD, CylinderU, CylinderV, Ei, EllipticCE, EllipticCK, 

  EllipticE, EllipticK, FresnelC, FresnelS, Fresnelf, Fresnelg, 

  GAMMA, GaussAGM, GegenbauerC, HankelH1, HankelH2, HermiteH, 

  JacobiP, KelvinBei, KelvinBer, KelvinHei, KelvinHer, KelvinKei, 

  KelvinKer, KummerM, KummerU, LaguerreL, LegendreP, LegendreQ, 

  LerchPhi, Li, LommelS1, LommelS2, MeijerG, Shi, Si, SphericalY, 

  Ssi, StruveH, StruveL, WeberE, WhittakerM, WhittakerW, arccos, 

  arccosh, arccot, arccoth, arccsc, arccsch, arcsec, arcsech, 

  arcsin, arcsinh, arctan, arctanh, cos, cosh, cot, coth, csc, 

  csch, dawson, dilog, erf, erfc, erfi, exp, hypergeom, ln, 

  polylog, sec, sech, sin, sinh, tan, tanh, JacobiAM, JacobiCD, 

  JacobiCN, JacobiCS, JacobiDC, JacobiDN, JacobiDS, JacobiNC, 

  JacobiND, JacobiNS, JacobiSC, JacobiSD, JacobiSN, JacobiTheta1, 

  JacobiTheta2, JacobiTheta3, JacobiTheta4, JacobiZeta, 

  InverseJacobiAM, InverseJacobiCD, InverseJacobiCN, 

  InverseJacobiCS, InverseJacobiDC, InverseJacobiDN, 

  InverseJacobiDS, InverseJacobiNC, InverseJacobiND, 

  InverseJacobiNS, InverseJacobiSC, InverseJacobiSD, 

  InverseJacobiSN, JacobiAM, JacobiCD, JacobiCN, JacobiCS, 

  JacobiDC, JacobiDN, JacobiDS, JacobiNC, JacobiND, JacobiNS, 

  JacobiSC, JacobiSD, JacobiSN, JacobiTheta1, JacobiTheta2, 

  JacobiTheta3, JacobiTheta4, JacobiZeta, WeierstrassP, 

  WeierstrassPPrime, WeierstrassP, WeierstrassPPrime, 

  WeierstrassSigma, WeierstrassZeta, Zeta, Im, Re, abs, argument, 

  conjugate, signum, Dirac, Heaviside, abs, max, min, piecewise, 

  signum, Dirac, Heaviside, LambertW, MeijerG, Stirling1, 

  Stirling2, Wrightomega, bernoulli, euler, unwindK, BellB, 

  CompleteBellB, IncompleteBellB, HeunB, HeunBPrime, HeunC, 

  HeunCPrime, HeunD, HeunDPrime, HeunG, HeunGPrime, HeunT, 

  HeunTPrime, MathieuA, MathieuB, MathieuC, MathieuCE, 

  MathieuCEPrime, MathieuCPrime, MathieuExponent, MathieuFloquet, 

  MathieuFloquetPrime, MathieuS, MathieuSE, MathieuSEPrime, 

  MathieuSPrime, cos, cosh, cot, coth, csc, csch, sec, sech, sin, 

  sinh, tan, tanh, arccos, arccosh, arccot, arccoth, arccsc, 

  arccsch, arcsec, arcsech, arcsin, arcsinh, arctan, arctanh, 

  fourier, fouriercos, fouriersin, hankel, hilbert, invfourier, 

  invhilbert, invlaplace, invmellin, laplace, mellin]

For more functions (such as "simplify", "factor", ... for instance), this could help you
?index/functions

Note there is also ?index/packages to get the list of all packages.

The problem is hard for two completely different reasons:

1. Your expressions are extremely complex.
2. And you have to defined what you consider as "common expressions".

Here is a slight variant of your toy problem

Note that Maple never"naturally" writes (d+h+k) nor 1*(d+h+k)
so your "y" appears as y = d+h+k+h*f+4*d 

x =  4*(2*d+h+k)+(d+h)*z+(d+k)*u
y =  a*(2*d+h+k)+(d+h)^2*q+(d+k)^2*p

What is better for you:
    (1) to say that there are 3 common sub-expressions 2*d+h+k, d+h, d+k and rewrite
        x = 4*A+B*z+C*u
        y = a*A+B^2*q+C^2*p
    (2) to say there are twos: d+h and d+k and rewrite
        x = 4*(B+C)+B*z+C*u
        y = a*(B+C)+B^2*q+C^2*p

Another example (e1, e2, ... are expressions of `+` type for instance
   x = (e1*e2 + e2*e3 * e3*e1) / (1+e1*e2*e3)
   y = (e1*e2^3 + e2*e3^2 * e3*e1^2) / (1+(e1*e2*e3)^2)
What would be for you the more concise expressions of x and y?

Let's go back to your problem now.

I only looked to the simpler "g1 vs g_1" case.
First thing: g_1 contains a new name "t" in a lot of expressions which, had this "t" name removed g1 and g_1 could appear quite close the one from the other; but it is no longer true of "t" remains in g_1.
So this "t" limits your expectations to concise expressions with several shared sub-expressions.

The first step (IMO) is to identify all those sub-expressions in g1 (let's say they are gathered into the set SE1) and in g_1 (the set SE_1) and build the intersection I of these two sets.
One issue is that I contains quite complex sub-expressions and very simple ones such as w^2 and w^3.
Replacing those latter doesn't seem to have no interest at all (what would be the point in writting w^2=A and w^3=A^(3/2)?), so it seems better to remove from I all expressions of type name^constant.

But even after doing so there remains expressions of large length (up to 120 or something alike) which indeed could deserved to be replaced by some new simple name, and much shorter ones like c+w.
Do you want them to be rewritten as, let's say A, or do you prefer to keep them unchanged?

The common expressions which deserve to be replaced by some simple name seem to me to be those whose length is much larger than the length of the "alias" name.

Here is an example An_attempt_mmcdara.mw

Personally I'm not in favor of extremely complex and unreadable expressions whose, IMO, the sole interest is to impress.
So I indulged myself in giving you some tricks to have more concise expressions.

Find-U-in-PDE_mmcdara_1.mw

One more time it is impossible to upload a file!

Waiting for this issue to be fixed, here are the key elements.

# The message is clear: plot can't be done sor the function to plot likely contains more
# indeterminates than x and y.
#
# Let's check that:

indets(eval(U(x, y, a, b, alpha, beta, lambda__1, lambda__2), t = -50))

                       {x, y, w[1], w[2]}

# So if you want to plot U, give w[1] and w[2] numeric values

More of this, what is yiur purpose writting

U :=(x, y, a, b, alpha, beta, lambda__1, lambda__2) -> rhs(eqt1)

while rhs(eqt1 doesn't depend on a, b, alpha, beta, lambda__1, lambda__2 ????
At the evidence you copied-pasted a bunch of commands from an answer you got (likely janhardo) without thinking.

If you want to define U correctly write this

# checking
indets(rhs(eqt1), name)
                     {t, x, y, w[1], w[2]}

# defining
U :=(x, y, t, w[1], w[2]) -> rhs(eqt1)

This will make explicit that contoutplot(U :=(x, y, t, w[1], w[2]), x=..., y=...) can't work as t, w[1], w[2] are still undefined.

I don't know if you are familiar with sensitivity analysis but the expressions you use are debatable.

Maybe they are of common use in epidemiology?

Download sensitivity_mmcdara.mw

discont seems to be a dead end, but fdiscont works well:

fd := fdiscont(ux, x=-10..10):

hranges := [ -10..op(1, fd[1]), seq(op(2, fd[i])..op(1, fd[i+1]), i=1..numelems(fd)-1), op(2, fd[-1])..10 ]:
sranges := map(r -> map(SFloat, r), hranges):
print~(sranges):
                  -10. .. 1.34837367547278353
           1.34924947400923112 .. 2.01735871146306556
           2.01825579287648305 .. 3.54283534891337881
           3.54368782900890489 .. 7.43153772737138940
                   7.43254178967007118 .. 10.

plots:-display( seq(plot(ux, x=r), r in sranges) )

singular-interval_mmcdara.mw