To get a Crout decomposition, first get a standard LUDecomposition of the transpose of the coefficient matrix, then transpose each of the returned matrices, then reverse the order that they are multiplied. Although that may seem like a lot of steps, it can all be done in one short command in Maple:

Download Crout.mw

Note the order specified in the output option: ['P', 'U', 'L'] rather than the usual ['P', 'L', 'U'].

You are trying to solve 2 equations for 1 variable. It's not usually possible to get an unconditional solution when there are more equations than variables being solved for. Here are 2 options:

• Add a 2nd variable (of your choice) to solve for,
• Accept a conditional solution.

To add a second variable (for example, Clm):

solve({tau = `τopt`, lambda = `λopt`}, {w, Clm})

To get a possibly conditional solution, use eliminate instead of solve:

eliminate({tau = `τopt`, `λopt` = lambdaopt}, {w})

The return will be a list of two sets. The first set is the solution(s) for w (only one in this case). The second set is expressions (only one in this case) that must equal 0 in order for the first set to be valid.

You've made a quite interesting discovery (Vote up): Overloads in modules don't explicitly have option overload, although they do otherwise contain the correct programmatic structure of an overload. So, the workaround for your procedure Isee is to explicily insert option overload when that's needed, like this:

Isee:= proc(a::procedure)
local vp:= interface('verboseproc'= 3), p:= eval(a);
    if op([5,1], ToInert(eval(a)))::specfunc(_Inert_OVERLOADLIST) then
        p:= subsop(3= overload, eval(a))  #op(3,a) is the explicit options
    fi;
    printf("%P", eval(p));
    interface('verboseproc'= vp);
    return
end proc:

There are several things that I did to make your double and quadruple summations faster. The most important and the most obvious is that any factor in a sum that doesn't depend on the summation variable should be factored out. Otherwise, those binomials (and other things) get recomputed many times,. Some other tricks:

1. GAMMA(s1+h+1)/(GAMMA(s1+1)*h!) = binomial(s1+h, h), even when s1 isn't an integer.
2. All necessary binomials can be pre-computed once, via very simple recursive formulae:
   [1, seq[scan= `*`]((n+1)/k - 1, k= 1..n)] generates the sequence binomial(n,k) $ k= 0..n (i.e., the entire (n+1)^st row of Pascal's triangle). By changing the order of subtraction to 1 - (n+1)/k, we generate (-1)^k*binomial(n,k).
   [1, seq[scan= `*`](1 + s/k, k= 1..upto)] generates binomial(s+k, k) $ k= 0..upto.
3. (beta1 - 1)^h*(-1)^h = (1 - beta1)^h, and of course this is trivial to generate recursively.
4. The additions of fixed terms in the denominators can be pre-computed:
   e.g., DL:= x+y + b1*(s1+h) + b2*(s2+v).
5. The undocumented command inner(A, B) computes the inner (or "dot") product of two lists of equal length; i.e., inner(A,B) = add(A*~B) = sum(A[k]*B[k], k= 1..nops(A)).
6. The updating assignment operators +=, *=, and ++ work exactly the same as they do in the numerous other languages that have them:
   S+= a  <=>  S:= S+a
   B*= k  <=>  B:= B*k
   d++  <=>  (d+= 1) - 1, i.e., it increments d and returns the pre-incremented value,

Here is my worksheet with all my improvements and verifications that my methods give the same results as @dharr.

Download FasterQuadrupleSum.mw

You wrote:

• The combine function applies transformations which combine terms in sums,...

Unfortunately, that is just the most simplistic purpose of combine. In my mind, the primary purpose of combine is to make function arguments more complicated (such as changing func1(x) to func2(2*x)); expand does the reverse, changing 2*x to x. The 2 could be replaced by any integer.

This has always seemed intuitively clear to me, and I never use a menu to do something that can be done with a typed command. Perhaps those two things are related to each other.

You simply need to use Typesetting:-Typeset in the plot option:

labels= [y(x), Typesetting:-Typeset(diff(y(x),x))]

To get a decimal version of the real value of a, use

evalf(a[1])

instead of your eval(a[1]). The f in evalf stands for floating-point. A more straightforward alternative is to replace your solve command with fsolve:

fsolve(period_eqn, a)
              8058.994329

There is no problem with Maple itself here. The problem is only with MaplePrimes's ability to display the worksheet inline. Downloading the worksheet link produces exactly the file that you uploaded.

Here, I've downloaded, executed, and re-uploaded your worksheet. If you download this link, I think you'll see that it's equivalent to your original:

Download predator_prey.mw

First, I'll compactify the code that produces your static plot, and I'll add hooks for controlling the parameters with more finesse. In particular, note the parameters option to dsolve and the solitary call to numsols with the parameters option:

restart:
vars:= [w, x, y, z]:
sys:= {
    diff~(vars(t), t)[]=~ ((y*z)(t), a*w(t), (x*z + w)(t), (-x*y)(t)),
    (vars(0)=~ vars||~0)[]
}:
params:= [a, w0, x0, y0, z0]:
param_vals:= [2, -0.727367040, -0.728244724, -0.237753623, 0.014225402]:
numsol:= dsolve(sys, numeric, method= rkf45, 'parameters'= params):
col:= [red, magenta, cyan, blue]:
numsol('parameters'= param_vals):
plots:-odeplot(
    numsol, `[]`~(t, vars(t)), t= 0..3, color= col, legend= vars(t)
):

The next step is to abstract that into a procedure that produces that static plot for given parameter values:

OnePlot:= proc(param_vals, tf)
    numsol('parameters'= param_vals);
    plots:-odeplot(
        numsol, `[]`~(t, vars(t)), t= 0..tf, 'color'= col, 'legend'= vars(t)
    )
end proc
: 

And finally the Explore:

Explore(
    OnePlot(params, tf),
    'parameters'= [(params[], tf)=~ (1. .. 3., (-1. .. 0.) $ 3, 0. .. 1., 3. .. 9.)],
    'initialvalues'= [(params[], tf)=~ (param_vals[], 3)]
);

Like this:

seq(cat("%15s" $ k), k= 1..4);

To continue, you need to say what mathematically makes the difference between solid and dashed lines.

params:= {
    h= piecewise(
        z <= d+1,   1,
        z <= d+4,   1-delta/2*(1 + cos(2*Pi*(z - 1 - 1/2))),
        z <= d+6,   1 
    ),
    w0= -c*h^2/4 + 3/64*(b*c-4*a)*h^4 + 19/2304*b*(b-4*a)*h^6,
    w1= c/4 + (4*a-b*c)*h^2/16,
    w2= (4*(b*c-4*a)-b*h^2)/256,
    w3= b*(b-4*a)/2304,
    a= x4*S*Gr*sin(alpha)/(4*x1*x5),
    b= 1/Da + x3*M/(x1*(1+m^2)),
    c= Dp/x1,
    Dp= 96*x1/((6-b*h^2)*h^4)*(F + a*h/24 - 11/6144*b*(b-4*a)*h^8),
    x1= ((1-phi1)*(1-phi2))^(-5/2),
    x2= (1-phi2)*(1 - phi1 + phi1*Rs1/Rf) + phi2*Rs2/Rf,
    x3= shnf/sf,
    x4= (1-phi2)*(1 - phi1 + phi1*RBs1/RBf) + phi2*RBs2/RBf,
    x5= khnf/kf,
    shnf= sbf*(ss2+2*sbf-2*phi2*(sbf-ss2))/(ss2+2*sbf+phi2*(sbf-ss2)),
    sbf= sf*(ss1+2*sf-2*phi1*(sf-ss1))/(ss1+2*sf+phi1*(sf-ss1)),
    khnf= kbf*(ks2+2*kbf-2*phi2*(kbf-ks2))/(ks2+2*kbf+phi2*(kbf-ks2)),
    kbf= kf*(ks1+2*kf-2*phi1*(kf-ks1))/(ks1+2*kf+phi1*(kf-ks1)),
    RBs1= 8933*16.7e6,
    RBf= 1063*1.8e6,
    RBs2= 6320*18e6,
    kf= 0.492,
    sbf= 6.67e-1, ss2= 2.7e-8,
    sf= 6.67e-1, ss1= 59.6e6,
    ks2= 76.5, kf= 0.492, ks1= 401,
    phi1= 0.01, phi2= 0.02, alpha= Pi/4, m= 0.5, Da= 0.1, Gr= 5, 
    delta= 1, S= 0.5,  d= 1,
    F= 1, z= 1
}:
                                      
W1:= eval[recurse](w0+w1*r^2+w2*r^4+w3*r^6, params):

plot(
    eval~(W1, M=~ [2, 5, 7]), r= 0..1,
    color= [red, blue, green], axes= boxed,
    legend= typeset~(M^2=~ [2,5,7]),
    labels= [r, w(r,1)], labeldirections= [horizontal, vertical],
    axesfont= [helvetica, 16], labelfont= [helvetica, 16],
    caption= typeset(
        ``(seq(v = eval[recurse](v, params), v= (m, Da, Gr, S, alpha, phi1, phi2)))
    )
);

If I translate your f into 1D-input, it becomes this:

f:= (u[1], u[2])-> u[1]^2*diff(u[2], x);

It's simply a nuance of Maple syntax that indexed variables such as u[1] are not allowed to be formal parameters of procedures. There is a simple fix that allows the variables to still display as subscripted:

f:= (u__1, u__2)-> u__1^2*diff(u__2, x);

c__1 is an alias. Aliases are very different from ordinary variables (and I think it was a bad design choice to use them for dsolve). One of the weird things about aliases is that if they are used in a procedure, they must be defined before the procedure is defined, whereas an ordinary global variable only needs to be defined before the procedure is called. To verify this, re-enter the definition of foo (without restart), and re-call foo (shown below). Also shown below is a safer way to code foo by being agnostic about what the constant will be.

Your Example 2 Question: Solving for _C1 works, even though the ode has c__1 , why?

Because _C1 is an ordinary global variable, not an alias. c__1 is an alias for _C1. To see this:

_C1;
                       c__1

Your Example 3 Question: Forcing arbitraryconstants = subscripted it still does not work inside proc. Why?

c__1 is still an alias. The weird interaction between aliases and procedures supercedes everything else.

Download AliasedDsolveConstants.mw

When applied to a first-order nonlinear ODE, the term homogeneous has a completely different and standard definition than the one used for linear ODEs and referenced by @C_R . There is no connection between the two definitions. A quick Google search turns up many references to this definition of homogeneous for 1st-order nonlinear ODEs, and it's equivalent to the defintion given on the Maple help page that you referenced. So, I think that that quote from Advanced Engineering Mathematics with Maple is not telling the whole story.

Each of your three algebraic solutions of your ode for diff(y(x), x) can be put in the form F(y(x)/x) (using three slightly different Fs). To see this, after your PDEtools:-Solve command do

simplify(eval(rhs~([sol]), y(x)= x*v), symbolic);

You'll see that none of the three depend on x. (Note that v = y(x)/x.)

What benefit do you think you're getting by using alias? Aliases are quite weird, and if you're thinking of XX as just another variable, you'll likely make mistakes. Why not instead use

export XX:= Sqr1;

or

macro(XX= Sqr1);

?

It seems to me as if you just want an abbreviation for a longer command name. The export above will do that. The only possible reason I know of to use an alias is that it's an abbreviation that gets back-substituted on output.