@zishihuandi 

It's not really a system of PDEs. It's a single, simple one-dimensional heat equation with some complicated boundary conditions (complicated because they involve integrals of the unknown function). I don't know how to proceed from there, but this observation substantially reduces the apparent complexity from the original nine equations.

Can you show an example that doesn't use a restart between the withs and doesn't work correctly?

Your eq1 isn't an equation. Did you mean for it to be equated to 0?

@tomleslie I think that it's only easy if the color scheme is zgradient, as it is in your example. To do a custom colorscheme gradient based on another function requires something like this:

FuncGradient:= proc(P::specfunc([realcons,realcons], POINTS), f, minMAX::range(realcons))
local m:= lhs(minMAX), R:= rhs(minMAX) - m, p;
     COLOR(HUE, hfarray([seq((f(p[]) - m)/R, p= P)]))
end proc:

To apply this to the example that I used in my Answer, do

XY:= RandomTools:-Generate(list(list(float(range= -9..9), 2), 2^12)):
plots:-pointplot(XY, colorscheme= [FuncGradient, (x,y)-> sin(x), -1..1]);

(The plot that this produces is exactly like the one in my Answer.)

This does run faster than the code that I gave in my Answer.

For the OP to apply this to their situation, use

plots:-pointplot(XY, colorscheme= [FuncGradient, f, 0.1..1]);

What do you mean by the "roots" of a multivariate polynomial?

What does "run a two for" mean?

I'm probably not the typical Maple user, but, FWIW, I find the syntax :-keyword= keyword completely intuitive, and I find the syntax keyword= keyword counterintuitive. Indeed, the first time that I encountered a situation similar to yours (it was a recursive procedure with a keyword parameter), 
:-keyword= keyword was the first thing that I tried. 

If you read through Maple library code of recent years, you'll see that it's rife with examples of 
':-keyword'= value, even in cases where that colon-dash makes no difference.

@Doug Meade 

Why not the following single command?

ans:= GenerateMatrix(eqns, [x,y](t))[2];

The OP expressed concern about memory usage. Both of your two-command solutions prevent the garbage collection of the coefficient matrix.

Maple can't even differentiate HankelH1 with respect to its first argument, let alone get a series. Note that diff(HankelH1(x,y), x) returns unevaluated whereas diff(HankelH1(x,y), y) does not.

@dharr I just think of indets as a command that returns all substructures of a given type. As such, it is extremely useful. I always use a second argument. I don't think of it as having anything to do with indeterminates. The one-argument usage should be deprecated.

Note that you can solve an equation for exp(a*d) and you can assign a value to exp(a*d).

@dharr This is a nuisance, but constants such as Pi are also considered names. So, keeping with the OP's use of type assignable, it's safer to use

indets(eq[1], assignable(name));

or, doing all equations at once,

indets({entries(eq, nolist)}, assignable(name));

Any multigraph can be converted into an essentially equivalent simple graph, suitable for use with Maple's GraphTheory package, by adding vertices of degree two. The process is as follows: If two distinct vertices are connected by n edges, then place a new vertex in the middle of n-1 of those edges. If a vertex is connected to itself by n edges, place a new vertex in the middle of all n of those edges.

@Preben Alsholm 

Vote up. Your commands can be combined with zip, which is, I think, a somewhat unusual usage of zip:

plots:-spacecurve(
     zip(solve, [(x-2)/3 , (y-1)/4 , (z-3)/3] =~ t, [x,y,z]),
     t= -5..5, thickness= 3, color= red, axes= normal
);

I don't know the answer to your Question. This Reply is just to bemoan the fact that what you want was extremely easy in Classic. You just typed Ctrl-R and typed your math expression as if it were a Maple command. The Ctrl-R still works in Standard, but it seems to only allow trivial math expressions, not stuff like summations and integrals.

This is the Classic feature that I miss the most :-( 

@tomleslie I believe that what the OP wants is to have the plots take up more space on the screen---to actually be larger.

@Dina 

If you need to retain the original matrix, as in your situation, then in-place operation is not appropriate.

It seems like you're trying to implement Cramer's rule. Here's a procedure that I wrote for it:

CramersRule:= proc(
     A::Matrix(square),
     b::depends(
          And(
               Vector[column],
               satisfies(b-> if op(1,b)=op([1,1],A) then true else error "Size mismatch" fi)
          )
     )
)
uses LA= LinearAlgebra;
     try
          <seq(LA:-Determinant(<A[.., ..k-1] | b | A[.., k+1..]>), k= 1..op(1,b))> 
          /~ LA:-Determinant(A)
     catch "numeric exception: division by zero":
          print("No unique solution.")
     end try
end proc:

The expression which substitutes b for the kth column of A is <A[.., ..k-1] | b | A[.., k+1..]>. This is valid even for the first and last columns. This generalizes Kitonum's three examples into a single command.