The command you need is convert(f(x), fullparfrac, x), possibly with added options. See the help page ?convert,fullparfrac.

Maple has a "sum over all roots of a polynomial" feature that usually makes seeing the individual factors unnecessary. However, if that's not useful for you, try including the factor option.

You get much better simplifications by calling dsolve for each numeric instantiation rather than getting a generic solution and then substituting parameter values,

Download DataFRame.mw

If you want to use Dimension, you just need to change uses RationalTriigonometry to uses LinearAlgebra. If you want to use type, you need to put 'Matrix' like that, with unevaluation quotes. That's because Matrix(1,3) is an executable command, not just a type.

(f,g):= <<a|b|c>>, <<d,e,h>>:
type(f, 'Matrix'(1,3)), type(g, 'Matrix'(1,3)), type(g, 'Matrix'(3,1));
                       true, false, true

Unevaluation quotes are not needed for types used in a procedure header:

f:= proc(p::Matrix(1,3)) ... end proc:

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 = `&tau;opt`, lambda = `&lambda;opt`}, {w, Clm})

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

eliminate({tau = `&tau;opt`, `&lambda;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);