Representing  (1 - 1e-18) requires 18 significant digits, whereas representing 1e-18 requires only 1 significant digit. So I find the following way much easier than all the others presented in this thread, requiring no increase in Digits from its default. I apply the following obvious identity:

Quantile(Normal(mu,sigma), 1-x) = mu - sigma*Quantile(Normal(0,1), x).

Thus

Quantile(Normal(0,1), 1 - 1e-18) = -Statistics:-Quantile(Normal(0,1), 1e-18);
      8.75729034878321

The following procedure will apply the rule only to those lns that occur as a result of the integration:

REALINT:= proc(f::algebraic)
local OrigLn:= indets(f, specfunc(ln));
   evalindets(
      int(args), 
      specfunc(ln), 
      L-> `if`(L in OrigLn, L, ln(abs(op(L))))
   )
end proc:

Example:
REALINT(1/(x*ln(x)), x);

If you're willing to have Maple tacitly make the necessary assumptions that'll get you that most simplification, you can use simplify(..., symbolic). Although I've never seen this occur in a practical situation, it is possible to construct expressions where the necessary assumptions are contradictory, and thus the simplification is not valid for any possible valuations.  See ?simplify,symbolic.

Manipulating square roots is difficult. Here's one tip: Almost all sqrt(...) are converted to (...)^(1/2), so searching expressions for sqrt won't work.

Here's a procedure to integrate recursively until either there's a term with no derivative or we get stuck with an unevaluated integral:

MyOdeInt:= proc(ODE::algebraic, Z::name:= NULL)
local 
   ode:= frontend(expand, [ODE]), #Distribute, but don't expand functions.
   #Find the integration variable:
   z:= `if`(Z=NULL, map2(op, -1, indets(ode, 'specfunc'(diff)))[], Z),
   r 
;
    if ode::`+` then #There's more than 1 term.
       if remove(hasfun, ode, diff) <> 0 then return ode fi 
    #single-term case:
    elif not hasfun(ode, diff) then return ode 
    fi;
    r:= int(ode, z);
    `if`(hasfun(r, int), r, thisproc(r))
end proc:

And you'd use it:

MyOdeInt(odetemp);

I tested on higher orders:

MyOdeInt(diff(odetemp, z$9));

which returned the same thing as the first command.

At the beginning of your worksheet, do

with(Units); with(Standard);
UseUnit('J/(mol*K)'):
UseUnit('J/(kg*K)'):

Then computations whose results can be expressed in those units will be. Computations include arithmetic, evalf, and combine(..., units). The raw extraction of R from ScientificConstants will stay the same, but you can force it to the preferred units with simply combine(..., units) rather than needing convert.

HFloat stands for "hardware float". These are the 64-bit representations of floating-point (or decimal) numbers that are nearly universally used by computer chips. (They're called "double precision" in many languages.) Included are representations of infinity, -infinity, and undefined.

For most purposes, you can think of HFloat(infinity) as simply infinity. If it appears as a result of integration, it indicates some type of divergence issue. Perhaps the integral is truly divergent. If it appears multiplied by some coefficient expression C, and C could possibly be 0, then you may have an indeterminate form, which'll require deeper analysis. There's no simple way to deal with this. It requires case-by-case analysis.

The results from solve should be put into a set. Then, iterate over the members (my s) of the set without regard for the index numbers (your k). Like this:

f:= x^2+2*x*y+2*y^3;
V:= [x,y];
gr:= VectorCalculus:-Gradient(f,V);
soln:= {solve({seq(gr)}, {V[]})};
(H,HD):= VectorCalculus:-Hessian(f, V, 'determinant');
seq(eval([H,HD], s), s= soln); #Change "s=" to "s in" if you like. 
#or
map2(eval, [H,HD], soln);

The reason for your error is that nops can't be applied to a sequence (unless it's a trivial sequence with 1 member)^[1]:

S:= 3,5,7,9: #a sequence
nops(S); #error
nops([S]); #make S a list
      4

So, it'd be a rare situation that nops(op(X)) would work at all! It could only work if nops(X) = 1. Also, note that nops(X) is by definition (see ?nops) the number of entries of [op(X)].

^[1]: Indeed, there's only a handful of cases where a nontrivial sequence can be used as a single argument to a procedure, and those are all built-in procedures.

It seems to be a bug. Here's a workaround:

J1:= int(exp(-t)/(1-t), t= 0..1-epsilon) assuming epsilon > 0, epsilon < 1;
J2:= int(exp(-t)/(1-t), t= 1+epsilon..infinity) assuming epsilon > 0;
simplify(J1+J2);
limit(simplify(J1+J2), epsilon= 0, right);

Yielding Ei(1)/exp(1), a real value. If the limit is incorrectly taken from the left, a related nonreal value is returned. Perhaps this has something to do with the bug.

Here's an implementation of procedure evalpw, as mentioned by Joe. I believe that I've covered every possibility; if not, let me know.

evalpw:= proc(PW, V::{`=`, {set,list}(`=`)})
description 
   "reduces a piecewise expression by considering its conditions" 
   " under evaluation without evaluating the branches"
;
option `Author: Carl J Love <carl.j.love@gmail.com> 3-Aug-2018`;
(* 
The 1st arg can be anything. If it's not a piecewise,
then it's returned unchanged. The 2nd arg is exactly like 
eval's 2nd arg, but it's only used to temporarily evaluate 
the conditions.  

The conditions are evaluated *in order* using PiecewiseTools:-Is
(undocumented, but its code is easy to read). If any evaluates
true, its branch is immediately returned. Any that evaluate false are 
removed from further consideration along with their branches. If all
evaluate false, then a default or "otherwise" clause is returned.
Otherwise, a reduced piecewise expression is returned.
*)
local  
   pw:= [op(PW)], n:= nops(pw), 
   k, #index to operands of PW  
   b, #boolean evaluation of condition
   R:= table() #undecidable conditions and their branches
;
   if not PW::'specfunc'(piecewise) then return PW fi;

   #Iterate over the conditions:
   for k by 2 to n-1 do
      b:= PiecewiseTools:-Is(eval(pw[k], V));
      if b = true then return pw[k+1] fi;
      if b <> false then (R[k],R[k+1]):= (b,pw[k+1]) fi;
   od;

   if numelems(R) = 0 then #No true, all false
      #If n::odd, use "otherwise", else use default.
      #Note that the default may be an index on 'piecewise' 
      #(it's an undocumented feature).
      return `if`(n::'odd', pw[n], op(0,PW)(false, ``))
   fi; 
   
   #Build residual piecewise:
   if n::'odd' then R[n]:= pw[n] fi;
   op(0,PW)(entries(R, 'indexorder', 'nolist'))         
end proc:

Make the first line of code of your worksheet:

restart;

This line should be in its own execution group.

If that doesn't solve your problem, then you'll need to attach your worksheet to a post, which you should usually do anyway when asking a Question. Use the green up-arrow on the editor's toolbar.

Elementwise operators (those that have ~ appended to their name) behave differently when both of the operands are container objects (such as lists). Observe:

f~([a,b], [c,d]);
       [f(a, c), f(b, d)]

Whereas you were hoping for [f(a, [c,d]), f(b, [c,d])]. To achieve that, you need to use map:

map(f, [a,b], [c,d]);

or, in your case,

map(coeffs, %, [x[1], x[2]]);

The above principle also needs to be applied to the preceeding collect command.

There's possibly an additional complication in this case: coeffs generally returns a sequence as output rather than a list, and you likely need some way to distinguish which coefficients came from which polynomial. In that case, use

map([coeffs], %, [x[1], x[2]]);

Then coefficients of each polynomial will be put in their own list.

Be careful with using the word Vector (with a capital V) because it has special meaning in Maple and may confuse your readers. A Maple Vector is a container object very similar to a list. Your polynomials are in a list, not a Vector. It turns out in this case that had they been in a Vector, the above code would work the same, except that the output would be a Vector.

In 2D input, you need a space between if and the following left parenthesis. This is not needed in 1D input.

Also, the parentheses are not necessary in either form of input.

You can only get a trigonometric solution with real coefficients when the non-real roots of the monic (i.e, M=1) characteristic equation occur as complex-conjugate pairs. That can only happen when C and K (or C/M and K/M) are real.

Perhaps you're interested in something like this:

restart:
f1:= 3.579167 + 3.537500*x1 - 2.645833*x2 - 0.250000*x1*x1 - 0.012500*x2*x1 + 0.175000*x2*x2:
plot(
   [seq(eval(f1, x1= k), k= 3..7)], x2= 3..7, 
   labels= [grape, preference], labeldirections= [horizontal, vertical],
   legend= (ginseng=~ [$3..7])
);

The following at least shows a correct way to input these equations and your known solutions, and plot those solutions using a truncation of the series. I couldn't figure out the pattern of the coefficients for the first solution, so I hard-coded the first 4 coefficients.
 

Download kdv_plot3ds.mw