I converted your Post to a Question.

It's extremely difficult to for us to re-execute your code in its current format. It's also difficult to follow references to equation numbers, such as the all-important equation 19. Please upload the worksheet using the green up arrow on the toolbar.

There's a mysterious symbol uminus0 or `&uminus0;` that begins to appear in your output as of MatrixInverse(XprimedX). I think that this has something to do with your trouble. Symbols with names like that are special mathematical symbols that come from the palettes and are often inadvertently pasted into code. There's no way to tell whether this is actually the cause of your trouble without close inspection of the actual worksheet.

Great Question; vote up. I don't have an Answer, but this is something that I've thought about before: How can I create my own parametrized custom distribution MyDist such that it can be invoked MyDist(a,b) in any situation where one could invoke, say, Normal(mu, sigma)? That is your goal, correct? If so, I think that you're well on your way to a method. I think that the key is constructing the module that is the 3rd attribute (see setattribute at ?attributes).

Here are five very small steps to help you on your way:

1. I see that you're using with, apparently simply to list the exports of the 3rd attribute. It seems to work, but I think that it's risky because with modifies the global state. You could replace with(A) with exports(A).

2. In all cases where you've defined a PDF, you've used both z and t, but you need to reduce that to a single name.

3. You've used ParentNames where you should've used ParentName.

4. No doubt you've read the examples at ?Statistics,Distribution, but they're worth re-reading carefully a few more times.

5. Avoid compound conditions in piecewise whenever possible: piecewise(a <= t and t < b, 1/(b-a), 0) can be and should be replaced by piecewise(t < a, 0, t < b, 1/(b-a), 0). I find that the latter form causes less problems, for example with symbolic integration. The only benefit that I know of the former form is that its prettyprinted display is more akin to standard mathematical typesetting. Maple's definition of piecewise guarantees that the leftmost condition that evaluates to true is the one that's used, which is what makes the two forms mathematically equivalent (for all real t). There's no need to construct a piecewise such that only one condition can be true because it's only the leftmost one that's true that matters.

@Christian Wolinski You're right: My g doesn't even do what I intended, which was that the output be (a,b)-> x-> a*x+b, with the internal a, b, and x being parameters, not globals. My testing was inadequate, and globals got through unnoticed.

You got a notice because you wrote an Answer (which used 'unapply' inside unapply), and I wrote a Reply to that Answer. At the very moment that I pressed Submit on that Reply, you deleted your Answer, but the notification had already been sent. The Reply itself was not posted because there was no longer any Answer to attach it to.

@tomleslie It's unfortunate that you suppressed the output of your most-important line. Had you not, it would've been obvious that your g re-executes unapply and eval on every invocation, which is undesirable.

@Christian Wolinski Although your command does ultimately produce a correct result (i.e., g(2,3)(4) yields 11), I don't think that it's what the OP really had in mind because the internal unapply is re-executed for every invocation of g.  

@ecterrab Thank you for pointing me to PDEtools:-dpolyform. It's an extremely powerful command. It's quite a shame that its help page---although extensive---doesn't show a single example where the command's input is anything other than algebraic---rather than differential---equations. Nonetheless, I was sucessful in applying it to all the ODEs in the OP's test suite.

However, dpolyform cannot be used to substantially reduce the complexity of my procedure ODEdegree (above) because, as you said, it returns an equation "that is polynomial in the unknown and its derivatives" [emphasis added] whereas the OP requires that only the derivatives, not the unknown, be put into polynomial form. Thus, usage of frontend is still needed to ignore the unknown but not its derivatives.

While I do not have expertise in the symbolic solution of differential equations, I fail to see the usefulness of the OP's concept of the "degree" of a differential equation. It does however seem to be a standard definition. My procedure implements that definition as per the OP's request.

@syhue Your "following steps" do not appear in your message. Either enter them into your message as plaintext or use the green up arrow on the toolbar to upload a worksheet.

If I may hazard a guess as to what you're working on: If you want to compute the remainder of (a^b)/p where a, b, and p are integers and b itself is very large, then do

a &^ b mod p

@666 jvbasha The command print can be used anywhere, in or out of a loop or procedure, to display output, including any type of plot.

@nm Thanks for the extensive testing. 

There's a command that can be used to test a procedure for a list of inputs and expected results: CodeTools:-Test. In this case you could do 

CodeTools:-Test~(
   ''ODEdegree''~((T1:= index~(tests, 1))), 
   index~(tests, 2), 
   label=~ convert~(T1, string)
);

Note the my procedure only works for syntactically correct ODEs, and it doesn't check that its input is such. Your 5th and 6th test cases are identical and have a syntax error: They have y where they should have y(x).

@Christian Wolinski 

I wasn't using the indets order; I was sorting by length. I think that I just need a tie-breaker for when the lengths are equal. For this, I'll choose set order. So, a small modification:

#Strict operand ordering:
S_O_O:= (x,y)-> length(x) < length(y) or length(x)=length(y) and x={x,y}[1]:

SubsIndets:= proc(e, t::type, T)
local J:= sort([indets(e, t)[]], S_O_O), r:= e;
   to nops(J) do (J,r):= subs(J[1]= T(J[1], _rest), [J[2..], r])[] od;
   r
end proc:
      
h := proc(x) global i; i := i + 1; T[i, op(0,x)](op(x)) end;
e:= f(f(f(f(x), g(y)), f(x)), g(y)):
i:= 0:  SubsIndets(e, function, h);
i:= 0:  subsindets(e, function, h);

Now the built-in version and my version return exactly the same thing (for your example, at least).

I agree. The left-side result seems "too recursive". I wrote my own versions of subsindets several different ways, and in every case, they retuned the same result for your left side and right side.

Please change the title of this Question. The banality and genericness of your title could lead to this Question being forgotten or overlooked. Something like "Serious discrepancy in subsindets/evalindets." (It affects evalindets also.)

@nm Here is a correction that uses frontend to temporarily ignore everything except the derivatives. I also added a check that the transformed ODE is a polynomial function of the derivatives. So, please test this on numerous examples:

ODEdegree:= (e::{algebraic, `=`(algebraic)})-> 
   try
      (e-> degree(
         e, 
         indets(
            e, 
            'typefunc'(
               (O-> `if`(O>1, 'typefunc'('`@@`'(identical(D),O)), specfunc(D)))
                  (PDEtools:-difforder(e))
            )
         )
      ))
         (
            ((e::satisfies(e-> degree(e, indets(e, specfunc(Diff)))::posint))-> 
               convert(e,D)
            ) 
               (frontend(
                  evala@Norm,
                  [convert(`if`(e::`=`, lhs-rhs, x->x)(e), Diff)], 
                  [
                     {
                         specfunc(Diff), `+`, `*`, 
                         'satisfies'(e-> hastype(e, specfunc(Diff)))^rational
                     }, 
                     {}
                  ]
               ))
         )
   catch:
      FAIL
   end try
;

You wrote:

• The hard part...is to see if the ODE can be written as a polynomial in derivatives....
  (Maybe Maple has some built-in functions to help with this algebra part?)

Better yet, Maple has a single command that does the whole polynomial conversion, evala(Norm(...)), so that part is not hard at all. It's in my code above, of course.

@nm Your indets[flat] seems okay, and better than my AllPowers procedure. I learned something new about the flat option (which is a fairly new option in Maple) from your example. I'm glad that you researched it! 

I hate repetition of code fragments, and I just figured out a way to avoid the repetition of identical(y):

indets['flat'](expr, identical(y)^~{1, algebraic});

Thank you. That's an impressive amount of formatting in a Maple document. Vote up. And I like how you show how Maple can be used for a large compendium of both non-mathematical and mathematical content. 

Some sugggestions:

1. I think that you should use standard left justification for the main exposition rather than the center justification that you used.

2. (I don't know if this was true in Maple 2015): Maple comes prepacked with all the geographical information needed to make a map of the World, which can be shown on a globe or projected to a flat map. You'd probably like the include something like that, probably in several places.

3. Surely Maple could be used to analyze remote sensing data (your 14th section). An example along those lines would be great. (My very first computer programming job (in 1979) involved analyzing satellite data of the ocean surface, primarily to detect chemical spills.) Maple comes prepackaged with all the "hooks" needed to download raw data from the Internet. 

@Carl Love Here is a procedure to address the modification needed based on your answer to Christian Wolinski's question:

AllPowers:= proc(e, y::algebraic)
local b,p;
   subs(
      {p= 1, b= y}, 
      indets(
         subsindets(subs(y= b^p, e), identical(b^p)^anything, e-> b^(p*op(2,e))), 
         identical(b)^anything
      )
   )
end proc:

Usage:
AllPowers~([y^2*sin(1/y) + y^(3/2) + x*y^7, y^2+y], y);

Note that the first argument can be any expression at all, including sets and lists, and the second argument can be any algebraic, not necessarily a name or function.