I don't know if that's the expected behavior, but if you add halt to the first event action and restart dsolve a few times with eventclear, you'll see that the solution moves forward only by 10^-14 or so before the event is triggered again. So one way to avoid this is to add discrete_variables=[d(t)] with the initial condition d(Tdec)=0 and change the first event to

[[x(t)=0, And(v(t)<0, t-d(t)>10^(-9))], [x(t)=0, v(t)=0, d(t)=t]]

Then the integration stops at CMAX.

Also, you can get rid of the range option. Again, I don't know if that's as intended, but it doesn't work with symbolic parameters. That is, dsolve isn't going to precompute the solution for the whole range when you set the parameters. The range option only messes up the behavior of eventclear and generates irrelevant warnings about extending the solution.

For prototyping, mod is quite convenient:

irr0 := proc(f, x, q) local a, n, d;
  n := degree(f, x);
  a := Powmod(x, q^n, f, x);
  if a mod q <> x then return false end if;
  for d in NumberTheory:-PrimeFactors(n) do
    a := Powmod(x, q^(n/d), f, x);
    if Gcd(a-x, f) mod q <> 1 then return false end if
  end do;
  true
end proc:

irr0(x^5+x^3+1, x, 2);
                              true

This looks almost the same as the pseudocode for the algorithm. There is some overhead, because the input has to be parsed every time. For the same reason I'm not sure how efficient Domains is. modp1 works directly with the lists of coefficients:

irr := proc(f) local modq, Equal, q, a, n, d, x, one;
  q := modp1(Modulus(f));
  modq := e -> modp1(e, q);
  Equal := (a, b) -> IsZero(Subtract(a, b));
  n := modq(Degree(f));
  x := Monomial(1, Indeterminate(f));
  one := One(Indeterminate(f));
  a := Powmod(x, q^n, f);
  if not modq(Equal(a, x)) then return false end if;
  for d in NumberTheory:-PrimeFactors(n) do
    a := Powmod(x, q^(n/d), f);
    if not modq(Equal(Gcd(Subtract(a, x), f), one)) then return false end if
  end do;
  true
end proc:

pp := modp1(ConvertIn(x^5+x^3+1, x), 2);
irr(pp);
                              true

Not as readable anymore, but note that you don't have to wrap every operation in modp1. Also, you cannot compare zppoly objects with =, your code fails because it returns false for any two different objects, even if they represent the same polynomial.

To convert a list of coefficients to zppoly, just use ConvertIn(list, x). There is no function to convert from zppoly to a list (I think there should be). A way that works but is not documented is simply op(zppoly). A way that's guaranteed to work is coeffs of ConvertOut.

comp := proc(g, h, f) local modq, q, A, B, BA, hi, ri, n, m, b, x, i;
  q := modp1(Modulus(f));
  modq := e -> modp1(e, q);
  n := modq(Degree(f));
  m := ceil(sqrt(n));
  B := Matrix(m, m, [ListTools:-LengthSplit]([op](g), m));
  A := Matrix(m, n);
  hi := Vector(m+1);
  x := Indeterminate(f);
  hi[1] := One(x);
  for i to m do
    A[i] := Vector([op](modq(Rem(hi[i], f))));
    hi[i+1] := modq(Multiply(hi[i], h))
  end do;
  BA := B . A;
  b := Zero(x);
  for i to m do
    ri := modq(ConvertIn(convert(BA[i], list), x));
    b := Add(b, Rem(Multiply(ri, Powmod(hi[m+1], i-1, f)), f))
  end do;
  modq(b)
end proc:

comp(pp, pp, pp);
                            1 mod 2

Look carefully at the outputs. PDE shows as PDE := M->{pde2, ...}. The variable pde2 doesn't have an assigned value. The reason is that after you press Enter on the pde1:=... line, the cursor jumps to the pde3:=... line. I have no idea why that happens. pde2:=... doesn't have syntax errors, and it has Executable Math checked. There seems to be some tiny invisible cell in front of it though. Just use Worksheet Mode, not Document Mode.

There also seems to be an error in the code. You have a value assigned to Pr, and then you use Pr as an index. Change Plots[Pr]:=... to Plots[M]:=... and Colours[Pr] to Colours[M]. Then the code works without errors, but I don't know if it works as intended.

I would say it's a bug. Integrating g1_inv as a black-box function works:

evalf(Int(v1 -> 1/(1+g1_inv(v1, .5)^2), 0 .. 1));
                          0.8903158471

If P(x, y) is the monic polynomial of degree 9, then this will be about three times faster than the same code using fsolve(P(x, y)) instead of Eigenvalues:

m := unapply(LinearAlgebra:-CompanionMatrix(P(x, y), lambda), [x, y]):
f := proc(x, y) option remember;
  LinearAlgebra:-Eigenvalues(m(x, y))
end proc:

plot3d([seq]((i -> (x, y) -> f(x, y)[i])(i), i = 1 .. 9), 0 .. 2*Pi, 0 .. 2*Pi);

In other words, if you want to compute the eigenvalues, just apply Eigenvalues to your original matrix.

You're giving fsolve a nonlinear system. Something like fsolve({x*2^GAMMA(a) = 4, a = 3}) already takes noticeable time. Instead, you should do eval(eqs, alpha=.4). Then, like you say, you have a linear system, which can be solved instantly. eqs can even be solved symbolically, with alpha being a parameter.

If you're curious about why this identity is true, it can be seen as follows. The Fourier transform of sinc is, up to a coefficient (u.t.c.), a unit pulse, or the pdf of the uniform distribution on -1..1. The transform of the nth power is  (u.t.c.) the convolution, or the pdf of the sum of n uniformly distributed variables.

To construct this pdf, take the uniform distribution on 0..1 instead. The Laplace (for a change) transform of the pdf is (1-exp(-p))/p. The pdf of the sum is the inverse transform of ((1-exp(-p))/p)^n, or the sum for j from 0 to n of

inttrans:-invlaplace(binomial(n, j)*(-1)^j*exp(-j*p)*p^(-n), p, t) assuming n::posint, j::nonnegint:
simplify(convert(%, factorial));
             -(-1)^(j+1)*n*Heaviside(t-j)*(t-j)^(-1+n)/(factorial(j)*factorial(n-j))

Getting back to sinc, its integral is the value of the Fourier transform at zero (u.t.c.), which is the midpoint of the support of the pdf, thus the same (u.t.c.) as the value of the above at n/2 (understanding the pdf as a continuous function), thus the sum from 0 to floor(n/2):

Pi*n*sum((-1)^j*(n/2-j)^(n-1)/(j!*(n-j)!), j=0..floor(n/2))

which is exactly your formula generalized to include odd n.

There is an example here, with a very simple mask.

It can be evaluated in terms of MeijerG. After a change of variables, the integrand is a sum of tau^sigma*MeijerG(a*tau)*MeijerG(b*tau) terms:

gr:=(beta,r)->(1/8*I)*(-HankelH1(1, beta*r)-(2*I)*BesselK(1, beta*r)/Pi)/beta^2;

subs(s = sqrt(tau), (1/2)*tau^sigma*convert(cos(alpha*s)*gr(beta, s), MeijerG, include = cos));   
           (1/16*I)*tau^sigma*sqrt(Pi)*MeijerG([[], []], [[0], [1/2]], (1/4)*alpha^2*tau)*(
           -MeijerG([[], []], [[1/2], [-1/2]], (1/4)*tau*beta^2)-
           I*MeijerG([[], [-1]], [[1/2, -1/2], [-1]], (1/4)*tau*beta^2)-
           I*MeijerG([[], []], [[1/2, -1/2], []], (1/4)*tau*beta^2)/Pi)/beta^2

and the result of integration of each one is a single MeijerG. But in this case the standard convolution formula gives

-(1/16*I)*sqrt(Pi)*((1/4)*alpha^2)^(-sigma-1)*
MeijerG([[-sigma], [-sigma-1/2]], [[1/2], [-1/2]], beta^2/alpha^2)/beta^2+
(1/16)*sqrt(Pi)*((1/4)*alpha^2)^(-sigma-1)*MeijerG([[-sigma], [-1, -sigma-1/2]], [[1/2, -1/2], [-1]],
beta^2/alpha^2)/beta^2+
(1/16)*((1/4)*alpha^2)^(-sigma-1)*MeijerG([[-sigma], [-sigma-1/2]], [[1/2, -1/2], []], beta^2/alpha^2)/
(sqrt(Pi)*beta^2)

and we want sigma=-1/2, for which two of the MeijerGs are not defined because the kernel then has a product of GAMMA(1/2+y) with GAMMA(-1/2-y). But, taking MeijerG([[-sigma], [-sigma-1/2]], [[1/2, -1/2], []], z), its kernel can be rewritten as

simplify(GAMMA(1/2-y)*GAMMA(-1/2-y)/(1/2+y))*GAMMA(3/2+y)/GAMMA(-y);
                       -GAMMA(-1/2-y)^2*GAMMA(3/2+y)/GAMMA(-y)

The expressions are identical, but now we get a fully functional MeijerG. Doing the same for the second problematic MeijerG, we get

ii := eval(-(1/16*I)*sqrt(Pi)*((1/4)*alpha^2)^(-sigma-1)*
MeijerG([[-sigma], [-sigma-1/2]], [[1/2], [-1/2]], beta^2/alpha^2)/beta^2-
(1/16)*sqrt(Pi)*((1/4)*alpha^2)^(-sigma-1)*MeijerG([[sigma], [-1, -sigma-1/2]], [[-1/2, -1/2], [-1]],
beta^2/alpha^2)/beta^2-
(1/16)*((1/4)*alpha^2)^(-sigma-1)*MeijerG([[sigma], [-sigma-1/2]], [[-1/2, -1/2], []], beta^2/alpha^2)/
(sqrt(Pi)*beta^2), sigma = -1/2);

The result checks out numerically:

evalf(subs(alpha = 1, beta = 3/2, [Int(cos(alpha*s)*gr(beta, s), s = 0 .. infinity), ii]));
              [-0.02356569227 - 0.03703703704 I, -0.02356569226 - 0.03703703703 I]

evalf(subs(alpha = 3/2, beta = 1, [Int(cos(alpha*s)*gr(beta, s), s = 0 .. infinity), ii]));
              [0.02364435567 + 0.04270509831 I, 0.02364435574 + 0.04270509831 I]

What does it simplify to?

Not sure how much can be done with that TIF image, I think that's just a magnification of what we see in the slide. https://i.imgur.com/7tXecNX.png is probably closer to the original image:

with(ImageTools): with(ArrayTools):

img := Read("https://i.imgur.com/7tXecNX.png");
h, w := op(convert(Size(img), list));

fft := SignalProcessing:-FFT(img); # in 2017.3, same as DiscreteTransforms:-FourierTransform

# make a better mask here
mask := Create(h, w, (y, x) -> piecewise(14000 < (x-w/2)^2+(y-h/2)^2 < 16000, 0, 1));

fft2 := Mask(CircularShift(Create(abs(fft), fit), floor(h/2), floor(w/2)), mask);
Embed(abs(fft2)^(1/4));

# restore abs(fft)*signum(fft), the imaginary parts should be small
fft2 := FitIntensity(CircularShift(fft2, ceil(h/2), ceil(w/2)), min(abs(fft)) .. max(abs(fft)));
img2 := Create(Re~(SignalProcessing:-InverseFFT(fft2*signum~(fft))), fit);

Embed([img, img2]);

Notice though that the original Matlab image is 256x256, as the dimensions must be powers of 2 for FFT :). If index.png is the image embedded in your post, then in Maple 2017.3:

with(ImageTools):

sz := 256:
img := Read("d:\\index.png")[140 .. 140+sz-1, 78 .. 78+sz-1, 1];

img2 := ln~(Array(abs(SignalProcessing:-FFT(img)), datatype = float[8], order = C_order));

img2 := FitIntensity(ArrayTools:-CircularShift(img2, sz/2, sz/2)):

PlotHistogram(img2, autorange);

Embed(Scale(img2, 2));

Transform img2 and then use CombineLayers to get the color function that you want.

One very unfortunate caveat is that if there is a mistake, such as forgotten ~ after ln, then Maple attempts to print the whole array in full in the error message, after which the only thing I can do is to close Maple from the Task Manager, losing all unsaved worksheets.

If sol is the output of simplify, then

applyrule(e::RootOf = 'applyop(simplify, 1, e)', rhs(sol));

There are simplify(...,symbolic) and combine(...,symbolic), which can do transformations that are not generically valid. Or you can pass assumptions to those functions to let the transformations go through. For positive values of the variables, the result of simplify(...,symbolic) is indeed a solution:

Download simp.mw

This is the numerator of the derivative:

rat := (1/4)*(-1/4+alpha*(-1+b)*e^2+((1-b)*alpha+(1/4)*b)*e)^2/((-1+e)^2*(b*e-1)*alpha*e*(-1+b));

ndrat := factor(numer(diff(rat, e)));
    ndrat := -(4*alpha*b*e^2-4*alpha*b*e-4*alpha*e^2+4*alpha*e+b*e-1)*(4*alpha*b*e^3+2*b^2*e^3-8*alpha*b*e^2-
      4*alpha*e^3+4*alpha*b*e+8*alpha*e^2-5*b*e^2-4*alpha*e+b*e+3*e-1)

No proof, but I suspect that:

either rat has two double roots on 0..1 and the maximum is between them, or there are no roots and no maximum on 0..1;

when there is a maximum, the corresponding root of ndrat is the second real root of the cubic factor (the first one will correspond to the maximum at negative e);

the transition happens when the two double roots merge into a single root of order 4.

You can use RootOf as a compact and numerically stable representation of the root.

e0 := RootOf(subs(e = _Z, op(3, ndrat)), index = real[2]);
    e0 := RootOf(-1+(4*alpha*b+2*b^2-4*alpha)*_Z^3+(-8*alpha*b+8*alpha-5*b)*_Z^2+(4*alpha*b-4*alpha+b+3)*_Z,
      index = real[2])

`αcrit`:= [solve](discrim(gcd(numer(rat), ndrat), e) = 0, alpha)[2];
    `αcrit` := (1/8)*(2*b-4-4*sqrt(1-b))/(-1+b)

Explore(plots:-display(
    plot(rat, e = 0 .. 1),
    `if`(`&alpha;crit` < alpha,
      plottools:-point([e0, subs(e = e0, rat)], symbol = solidcircle, symbolsize = 15),
      TEXT([.5, 1.5], "no maximum")),
    view = [0 .. 1, -.1 .. 2]),
  b = .1 .. .99, alpha = 1.01 .. 10);

Maple can solve some non-polynomial inequalities, so you can take the generated decomposition of the (x,y) region and feed it into int():

cad := solve({abs(abs(x-y)-abs(x+y)) >= 2*y-x+1, (x+1)^2+(y+1)^2 <= 2}, [x, y]);

f := (a,b) -> if a = y then ymax = b elif b = y then ymin = a elif a = x then xmax = b else xmin = a end if:

g := s -> if subs(s, {xmax, xmin, ymax, ymin})::freeof({xmax, xmin, ymax, ymin}) then
    int(subs(s, ymax-ymin), x = subs(s, xmin) .. subs(s, xmax)) else 0 end if:

(simplify@add@curry(map, g))(applyrule(
  ['`<`(a::anything, b::anything) = f(a, b)', '`<=`(a::anything, b::anything) = f(a, b)'],
  cad));
                -arcsin((3/13)*sqrt(2)+(2/13)*sqrt(11))-arcsin((4/17)*sqrt(15)+(1/17)*sqrt(2))+
                (1/663*(-51*sqrt(11)+39*sqrt(15)))*sqrt(2)+2*Pi+419/663