Two things needed: 

1. Replace U and V with U(x,t) and V(x,t) throughout pde1 and pde2.
2. The first argument to pdsolve should be {pde1, pde2}.

These are necessary to solve your problem but not necessarily sufficient. I can't get on my computer right now to check. (There's a "no laptops" sticker on the café table I'm at.)

How do you expect n! to be defined for negative integers n other than by extension to GAMMA?

I think that Maple is primarily designed for high and wide computational power rather than to be educational or didactic. That requires that most functions be defined over the largest part of C possible. But, you can restrict the domains using assuming. In this case, you should append assuming n::posint to the original product. The resulting formula will contain GAMMA functions, but it'll evaluate correctly for your range of n regardless of whether you convert to factorial.

It seems that your attempt to upload the worksheet didn't work. Could you try again please?

Are you trying x[1] or x__1 for the variable? If it's either of those, it's surely the first. And did you tell it the name stem that you wanted to use? In that case, tell it to use :-x (colon hyphen x). Or did you let it choose the name stem? Then it must be _x[1]. Unassigned global names that library procedures generate for output always begin with underscore.

x[1] is a name with an index subscript; x__1 is a name with an inert subscript. They display the same (as subscripts) in pretty-printed output.

If you have output that you don't understand how to type, often applying the command lprint to that output improves understanding:

lprint(%);

Use expand rather than simplify. That's just being specific about the type of simplification that you want.

Your sum stops at m. I think that should be n.

The best way to fix the mantissa is to set the Digits environment variable. Maple internal floating-point arithmetic is base 10 with an integer mantissa with maximum of Digits digits. If I wanted to show the same anomaly of floating-point computation, I'd use this loop:

HarmonicSum:= proc(
   digits::posint:= Digits,
   rounding::identical(0, nearest, -infinity, infinity, simple):= Rounding
)
local S0:= -1., S1:= -1., S2:= 0., k;
   UseHardwareFloats:= false; #Optional
   Digits:= digits;
   Rounding:= rounding;
   for k to infinity while S2 > S1 do
      (S0,S1,S2):= (S1, S2, S2 + 1./k)
   od;
   ['n' = k-2, 'S'[k-3]=S0 , 1/(k-2) = 1./(k-2), 'S'[k-2]=S1, 1/(k-1)=1./(k-1)]
end proc: 

I think that that makes it explicit what's going on, right? And it gives you control over the order that things are added, which is important for this example.

HarmonicSum(7);

The output is the number of terms, the penultimate partial sum, the final term, the final sum, and the first term not used because it didn't change the sum. Because of the base-10 arithmetic, it's really easy to see what's going on.

If x is a Maple float, then [op(x)] will return the mantissa and exponent as an ordered pair.

Update: I added an option to allow the user to set the rounding mode.

When you're composing a post on MaplePrimes, the toolbar in the editor has a green up-arrow. It's the third item from the right in the top row of the toolbar. Let me know if you have any further trouble with uploading.

This is all that you need to do: Somewhere after the end proc of the definition of procedure si but before your usage of si, do

si:= subs('s'= s, eval(si)):

That's all!

When subs is used on a procedure in this way, it has the profound, almost spooky or magical, effect of changing global variables into parameters or locals, whichever is appropriate for the situation. This usage of subs is not documented on its help page (?subs), and I don't know if it's officially documented anywhere, but it is extensively used in Maple library code, so I don't think that it's ever going away. This usage of subs provides an extremely powerful technique for dynamic code generation.

In actual usage, I'd keep si as a "procedure-schema" and use the subs to create a procedure of a different name, as in

si1:= subs('s'= s, eval(si));

That way you can re-use si if you change to a new s. Also, it's not necessary that the s inside the procedure and the s outside have the same name. Indeed, it'd preferable that they were different. I often see experts using a name like _DUMMY_s for the one inside the procedure.

This usage of subs is similar to using unapply. The limitation with unapply is that it only works on a single expression, not a whole procedure. From Acer's code and Joe's comment about it, you can see that it wasn't super-super-easy to mangle the if statement from your procedure into something that unapply would accept. (Well, I'm sure it was easy for them, but they're the two greatest Maple programmers in the world IMO.) If it had been a do loop, it may have been impossible. If you read the code of unapply (by doing showstat(unapply)), you'll see that it's using the subs-into-procedure trick anyway.

The Maple command is timelimit. See ?timelimit

Example:

MyComputation:= ()-> 
   isprime((3^nextprime(rand(2^20..2^30)())-1)/2)
:
try
   r:= timelimit(10, MyComputation())
catch "time expired":
   r:= "Time expired."; #And proceed without error.
catch: #Some non-time-related error occurred:
   error
   #The actual error message will be printed at execution time.
end try;

Applying the methods from either of the Answers above could produce a negative result if the data were different or if simply the order that they were presented was different. This is not an error, but you do need to take the absolute value. So, to be general, you should use

abs(TripleScalarProduct(P1-P5,P6-P5, P8-P5))

or

abs(V1 &x V2 . V3)

The triple scalar product can also be computed as a 3x3 determinant, so my symbol-laden entry is

P||(1,5,6,8):= 
   <-17/12,11/36,-65/36>, 
   <-11/12,-7/36,-47/36>,
   <-2/3,-1/9,-8/9>,
   <-5/4,-1/12,-17/12>
:
abs(LinearAlgebra:-Determinant(`<|>`(map(`-`, [P1,P6,P8], P5)[])));

Under the assumptions that I mentioned in a Reply above, the following works:

TypeTools:-AddType(
   hname,
   proc(n)
   local r:= n;
      if not n::name then return false end if;
      while not r::symbol do r:= op(0,r) end do;
      evalb(substring(r, -1..-1) = ':-h')
   end proc
):      
select(s-> andmap(hastype, remove(evalb, s), hname), Sa1);

Explanations:

The type hname returns true iff the item being type checked is a name whose "main stem" (or symbol) ends with h. A while loop is used to strip off an arbitrary number of index subscripts. This technique will not handle double-underscore subscripts because they do not negate the symbol type.

If s is a set or list of equations, then remove(evalb, s) removes equations of the form X=X. From what remains of s, we require that every equation contains at least one hname.

Yes, you can change file permissions with the system or ssystem commands, unless, I guess, there's some rule in force that prevents Maple or other arbitrary programs from changing permissions. In that case, you'd need to grant Maple permission to change file permissions.

Some of the other ideas presented above may be better solutions for your actual problem, but taking your question at face value, a direct solution is still possible. By "face value", I'm interpreting your question as

• How can I use symbolic constants in an expression so that they will remain symbolic at all times except in those situations, such as plotting, where numeric values are required?

To do this, for each constant C, define a procedure (that takes no arguments) named `evalf/constant/C` that returns a numeric value. (And note that Maple considers any explicit numeric constant as a procedure that returns itself.) So, for your example, you could do

`evalf/constant/F0`:= 1:
`evalf/constant/rho0`:= rand(0. .. 1.):

and then proceed to define your vector field and do the field plot exactly as you had been.

Note that the procedure for rho0 above returns a different random value each time that it's used. The procedures can be of arbitrary complexity, and can include selection from any set or statistical distribution, as long as they have no arguments and return explicit numeric values.

If you go for a random selection, I can't guarantee that multiple instances of the same constant in one expression will get the same value. I think that it's likely that they will (due to evalf's cache or remember table), but it's not guaranteed. Likewise, it's not guaranteed that separate evaluations will give distinct values; you may need to clear the cache to guarantee that.

Do

remove(s-> true in map(e-> (is(rhs(e) < 0) assuming positive), s), Sa1);

Note for the experts: ormap won't work above because it refuses to handle the FAILs returned by is.

Short explanation: Only the right side of each equation is examined because it's assumed that the left sides are all simple names, as would be the case if these sets were returned by solve. For each right side, the is command tries to prove (using the ordinary rules of deductive logic and elementary algebra that we all know) that the expression must be negative under the assumption that all its names represent positive reals. If that proof succeeds (as indicated by returning true) for any case, the corresponding set containing that equation is removed. As is well-known and oft-discussed, is's ability to prove such things is quite limited, and it often returns FAIL, indicating that it can neither prove nor disprove a proposition. If this paragraph indicates that I've assumed something that you didn't intend, please let me know.

Use mul instead of product and add instead of sum because these operators allow the indices k to come from sets instead of ranges. The index sets are easy to construct with $ and minus. Like this:

restart:
N:= 4:
A:= Matrix(
   N$2,
   (i,j)-> if i=j then
      add(1/(x[i] - x[k]), k= {$1..i-1, $i+1..N})
   else
      mul(x[i]-x[k], k= {$1..N} minus {i,j})/mul(x[j]-x[k], k= {$1..j-1, $j+1..N})
   end if
);

If you're using 2D Input, the if ... end if might need to be changed to `if`(...) like this:

A:= Matrix(
   N$2,
   (i,j)-> `if`(
      i=j,
      add(1/(x[i] - x[k]), k= {$1..i-1, $i+1..N}),
      mul(x[i]-x[k], k= {$1..N} minus {i,j})/mul(x[j]-x[k], k= {$1..j-1, $j+1..N})
   )
);

(I'm not sure because I never use 2D Input.) Regardless, this latter form will work in any of Maple's interfaces.

Here is a procedure for it.

GramSchmidt:= proc(B::{Vector, list}, IP)
local n:= numelems(B), R:= Vector(n), M:= Vector(n), i, j;
   for i to n do
      R[i]:= B[i] - add(IP(R[j],B[i],args[3..])/M[j]*R[j], j= 1..i-1);
      if R[i] = 0 then error "Linear dependence detected at vector", i end if;
      M[i]:= IP(R[i]$2,args[3..])
   end do;
   <R|M>
end proc:
      
GramSchmidt(
   [seq(x^k, k= 0..5)],
   ((f,g,x::name)-> int((1-x)^(3/2)*f*g, x= 0..1)),
   x
);

The output is in two columns. The first is your orthogonal polynomials. The second is the constants that you called h_n.