In the somewhat confusing terminology of combinatorics, a partition (of an integer, not a set) is just a composition without regard to order. So, the following works:

select(p-> nops(p)=3, combinat:-partition(7));

Furthermore, it is obvious that the maximum possible entry in any such partition is 5 = 7-3+1, so the above command could be made more efficient as

select(p-> nops(p)=3, combinat:-partition(7,5));

Inside your loop, you have the statement

p:= solve(egp, p);

The first time through the loop, it works. But after this statement is executed once, p is assigned a value, and hence it's no longer a symbolic variable. So an error occurs the second time the solve is executed.

To correct the problem. make the first statement inside the loop (after all the for lines):

p:= 'p';

That's how you remove any assigned value that a variable has. (And this statement safely does nothing the first time through the loop, when p doesn't have a value anyway.)

There may be some other problems with your code (the output looks weird), but the above issue is definitely the major issue.

Note that your loops as they're currently bounded will produce close to 100 million iterations (92,692,600 to be exact). Do you really want to print output for every iteration? 

If S is a set of sets, then the following all produce identical results:

1. `union`(S[]) (works for S::{list,set}(set))
2. op~(S) (works for S::set({list, set}))
3. {seq(s[], s= S)} (works for S::{list,set}(list,set))

All are quite fast, but I've often wondered which is fastest. I haven't been able prove it definitely. I've convinced myself that 3 is faster than 2, which is a disappointment since 3 uses annoying, stupid syntax. But I haven't been able to rank 1 in relation to 3. What these operations do is not entirely trivial, because Maple sets have duplicates removed and are stored sorted into a canonical order.

If the members of S are not themselves sets, then there's no replacement for `union`. For example, you can express the union of pure symbols: `union`(A,B,C).

I am curious what process leads to Maple returning an ODE solution of that form. But, taking your word for it, suppose that you have an expression (possibly an equation, as above, or any other expression) e and you want to extract from it all function names that begin with Z and are followed by an integer. Then do

map2(op, 0, indets(e, typefunc(suffixed('Z', integer))));

If the names actually begin with _Z, then simply change 'Z' to '_Z'.

The other Answers seem to be not seeing the forest because of the trees. They suffer from unneeded specificity and are a bit difficult to generalize. Suppose that you have a set (or list) of prefixes P and a set (or list) of suffixes S. Then what you want is a function (indexing or suffixing in this case) double-mapped over the setwise exponentiation S^P (i.e., the set of all functions (in tuple form) from P to S). It can be generated by

P:= [CPC, SIZE, SH]:  S:= [h, m, l]:
seq(rtable((1..nops(S)) $ nops(P), ()-> index~(P, S[[args]]), order= C_order));

The order= C_order only affects the overall order of the sequence of tuples and may be omitted if you don't care about that.

The above generates indexed names (i.e., of the form A[B]) rather than suffixed names (of the form A__B). If you'd rather have suffixed names, then change index^[*1] to ((p,s)-> p||__||s). But note that programmatically generated suffixed names are always global, no matter how they are created and regardless of the lexicality of their components.

^[*1]Indeed, index could be replaced by any function whose domain is P x S. The canonical choice would be `[]` (the sequence-to-list operator). Whatever function is used, it gets mapped over every tuple of S^P (which I referred to above as double-mapping over S^P itself). It's important to realize that the fundamental combinatorial structure being generated is S^P, and the double-mapped function is simply an afterthought.

Use 4 distinct exact numeric values of t (values of rand() should work well). Evaluate each of the Qs at a different t value, then use regular solve on the resulting system of 4 equations. If you have extraordinarily bad luck, there might not be a unique solution due to  duplication of an equation. In that case, just use another set of t values.

Another one-liner, completely different from John's, is

add~([entries(ListTools:-Classify(content, [op(p)]), nolist)]);

Suppose that you have a list or set of equations such as [x(t) = 5, y(t) = 10] and that list or set is assigned to E. Then all that you need to do is 

eval([x(t), y(t)], E),

and you get [5, 10].

Usually what is done is to form one equation for each coordinate dimension and then solve the resulting system of nonvector equations. If you elaborate or give a small example, I can say more.

Here's a module for simplifying any collection of any real intervals (expressed with RealRange and/or as discrete sets of points) combined arbitrarily with any of the set operators `union`, `intersect`, and/or `minus`.

SimplifyIntervals:= module()
local
   x::identical(x),
   ModuleLoad:= ()->
      TypeTools:-AddType( #recursive type: put base cases 1st!
         'RealIntervals', 
         Or(
            specfunc(RealRange),
            set(realcons), 
            'specop'('RealIntervals', {`intersect`, `union`, `minus`})
         )
      ),
   Deconstruct:= (IN::RealIntervals)->
      if IN::specfunc(RealRange) then convert(x::IN, relation)
      elif IN::set(realcons) then Or((x =~ IN)[])
      elif IN::`minus` then And(thisproc(op(1,IN)), Not(thisproc(op(2,IN))))
      else `if`(IN::`union`, Or, And)(map(thisproc, [op(IN)])[]) 
      fi,   
   ModuleApply:= proc(IN::RealIntervals)
   local pts, ints;
      (pts,ints):= selectremove(type, [solve(Deconstruct(IN), x)], realcons);
      `union`(sort(ints, key= (p-> `if`(op(1,p)::realcons, op(1,p), op([1,1],p))))[])
  	  union {pts[]}       
   end proc
;
   ModuleLoad()
end module:

Example:

IN:= (((RealRange(-infinity, 1) intersect RealRange(Open(-1), infinity))
   minus RealRange(-1/3, 1/3)) minus {0}) union {-2,2};

SimplifyIntervals(IN);

(I tried to make a procedure `simplify/intervals`, but simplify just ignored it when I used simplify(IN, intervals). Any idea why?)

You seem to have ignored the z column of your matrix A. So your 3 planes can't be correct if they're supposed to represent A.<x,y,z> = 0.

Here's a way that does the computations and plots the planes and the line of intersection all in one command:

A := <1, 1, -2 | 3, 4, -7| -5, -8, 13>;
plots:-display([
   plot3d(
      #Solve each row for z using matrix operations:
      [seq(-A[..,..2].<x,y>/~A[..,3])], x= -8..8, y= -20..20, 
      plotlist= true, color= [blue, red, green], style= surface, transparency= [.3,.3,.1]
   ),
   #The line:
   plots:-spacecurve(LinearAlgebra:-NullSpace(A)[]*t, t= -3..3, color= red, thickness= 3)
],
   axes= frame, tickmarks= [[0=0]$3], labels= [x,y,z]
);

 LinearAlgebra:-NullSpace is essentially the same as your LinearSolve command, but returns the result in a more convenient form (a basis).

You had the right idea; there are only some minor syntax corrections needed. Try this:

restart:

Eqset:= {diff(x(t),t)=(-x(t)+x(t)^2-3*x(t)*y(t)), diff(y(t),t)=(-y(t)-5*x(t)*y(t)+y(t)^2)};
#The above line is your code, unaltered.

V:= x^2/2 + y^2/2;  #I removed the braces { }.
dV:= expand(eval(diff(V(t),t), Eqset)); #derivative wrt t (b/c Maple knows the symbolic chain rules)

plot3d(dV,x=-0.00005..0.00005,y=-0.00005..0.00005); #Your command, unaltered, now works.

(Today I learned... that plot3d allows functional variables such as x(t) and y(t) to be treated as scalars x and y.  Many Maple commands allow this, but I wasn't aware of its use in plotting commands.)

As you've suggested, it's theoretically possible to find some inequality involving x and t that'll guarantee that the other part is at most 10. But, why bother, unless you have some other reason for finding that inequality? You can just use min to explicitly set the maximum at 10, like this:

Y:= (x,t)-> min(10, -(11*sin(t/2+1.145)-10)*((x-37)^2)+(-3.5*sin(t)+ 1));

I think that it's best to make it a function of both x and t. If you insist on it being a function of x alone, then just remove the t from the left side of the arrow. 

Note that it's not necessary (or desirable) to use unapply to create this function because there's no meaningful evaluation possible before values are given for x or t. But neither is it incorrect to use unapply. It's just a matter of using a tool that's much larger than what's required to do the job.

It should be done with fsolve rather than solve because fsolve gives you more control over the range. Like this:

restart;
Digits := 32;
t0 := 1;
eq := 1-w*v^2-2*v*exp(-t/v);
equ := eval(eq, v = -t/ln(u));
#Lines above are exactly your code.

us:= subs(_E= eval(equ, t= t0), w= _W, _W-> fsolve(_E, u= -infinity..1));
vs := w-> -t0/ln(us(w));
plot(vs, 0 .. 10, labels= [w, ``], view= 0..1);

Do

plot([rhs(f), Y, Y= 0..10], labels= [typeset(epsilon/`2`/Pi), Y]);

Using plot rather than implicitplot is usually preferable when either the equation can be solved for one of its variables or its two variables can be expressed as functions of a third variable.