Try something like this:

A:=<a,b,0; 0,0,c>;

consts:=[a=2,b=3,c=4]:
A1:=eval(A, consts):
A2:=eval(A, c=66):
A1,A2;

simplify(expand(%));
collect(%,indets(%,name));  #optional

Among the methods used by fsolve there are:

For equations: Newton, Secant,  Dichotomic, inverse  parabolic interpolation.
For systems:  Newton, methods based on approximating the Jacobian and partial  substitutions.
Ostrowski's method (enhanced Newton) is also used sometimes.
If a method does not work, another one could be invoked.

In order to find what is effectively used, set
infolevel[fsolve]:=1;  #or greater if you want more details

Note that it is not difficult to find examples where fsolve fails but solve (or even the user, by hand) solves it.
E.g.
fsolve({(2*x+y+1)*exp(-(x+y-1)^2),   (3*x+2*y-1)*exp(-(x-y-1)^2)});

eq := z^2+y^3+x^4 - ln(x+y+z);
z__xy := implicitdiff(eq, z, x, y);
z__yx := implicitdiff(eq, z, y, x);
'z__yx' - 'z__xy' = z__yx - z__xy ;

with(plots):with(plottools):
display(disk([0,0], 1, color=red,transparency=0.5),disk([1,0], 1, color=green),scaling=constrained);

Probably you want:

with(plots):
ball := proc (a, b) fieldplot([x-a, y-b], x = -2 .. 2, y = -2 .. 2) end proc;
ball(1,1);

Your version makes no sense because in the procedure ball, x and y are is evaluated (substituted) with 1 producing
fieldplot([1, 1], 1 = -1 .. 1, 1 = -1 .. 1); #nonsense

We may use Poisson' formula:

uu:=x^2/2 + 1/(2*Pi)*(1 - (x^2+y^2)) *
      Int( 1/( (x-cos(t))^2+(y-sin(t))^2  )*(-cos(t)^2/2), t=0..2*Pi):
simplify(value(uu)) assuming x<1/2,y<1/2,x>-1/2,y>-1/2;

Of course, the result is valid for all (x,y) not only for |x|<1/2, |y|<1/2 but unfortunately the assume facility is not able to do this.
Note that value(uu) without assumptions ==> undefined ...        [???]

Note also that pdsolve is not able to find the correct arbitrary functions _Fnn

You should use lists instead of sets, otherwise the order of the elements is unpredictable.

f:=[1,2,3,4]:
h:=[1,2,4,5]:
select(i->f[i]=h[i], [seq(1..4)]);

maybe:

f[%];

or

select(i->i=h[i], f);

You must have a mistake somewhere because your analytic solutions satisfy the system only when C=0 or J=0, i.e. the trivial case.

pdesys := {diff(C(x,t),t)-phi*epsilon*J(x,t)*C(x,t), diff(J(x,t),x)-epsilon*J(x,t)*C(x,t)}:
Jsol := (x,t)-> J__0*exp(J__0*epsilon*phi*t)/(exp(J__0*epsilon*phi*t)+exp(C__0*epsilon*phi*x)-1):
Csol := (x,t) -> C__0*exp(C__0^epsilon*phi*x)/(exp(J__0*epsilon*phi*t)+exp(C__0*epsilon*phi*x)-1):
subs({J=Jsol,C=Csol},  pdesys): simplify(%);

X := -0.6356560300e-1*ln(Y+200.+17.54410643*Y^(2/3)+102.5985568*Y^(1/3))/a+
         0.6356560300e-1*ln(Y+200.)/a-.2201977080*arctan(.1974510146*Y^(1/3)
         -.5773502693)/a+.1096187623/a:
vx := <1500, 1340.00, 1135.00, 982.00, 884.15, 704.72, 520.00, 287.00, 70.00, 0.>:
vy := <0., 4.28, 7.77, 7.30, 9.00, 13.00, 25.85, 28.91, 38.48, 50.00>:

Statistics:-Fit(X, vy, vx, Y,output=parametervalues);
    [a = 0.115891288738057e-3]

f:=Statistics:-Fit(X, vy, vx, Y):
p1:=plot(vx,vy,style=point):  p2:=plot( [f,Y,Y=0..50],color=red):
plots:-display(p1,p2);

In Maple 2016.1 (64 bit)  I obtain the correct answer:

[12, [x = 1, y = 2]]

This is normal because verify does a kind of syntactic check (extended by some options).

I am a bit surprised that the probabilistic check testeq also fails.

It would have been nice if j:=16117;N:=3;  were provided in the code.
Even nicer if the unknowns had reasonable names (for an outsider).
So, with a simple subs your system is:

sysX := [x6,
 x20-x23,
 x21-0.3505865589e-5,
 x22-0.5304364281e-5,
 x5-x14,
 x24*(-0.3915554290e-1*x12-0.1903748329e-1*x6+0.8260795999e-1)-3.876387504,
 x24*(-0.3860115660e-1*x12-0.1876793978e-1*x6+0.7836678184e-1)-2.040147478*10^6*x23,
 x24*(-0.1876794098e-1*x12-0.9892449327e-2*x6+0.3810204607e-1)-2.040147478*10^6*x9,
 x25*(-0.3915554290e-1*x13-0.1903748329e-1*x7+0.8260795999e-1)-3.876387504,
 x26*(-0.3915554290e-1*x14-0.1903748329e-1*x8+0.8260795999e-1)-3.876387504,
 x1-.9724029753*x6-x12,
 x10-.9724029753*x7-x13,
 x11-.9724029753*x8-x14,
 x3-x12+x10,
 x3-.25*x18-.25*x16,
 x4-x13+x11,
 x25*(-0.3860115660e-1*x13-0.1876793978e-1*x7+0.7836678184e-1)+2.040147478*10^6*x20-7.152482840,
 x25*(-0.1876794098e-1*x13-0.9892449327e-2*x7+0.3810204607e-1)+1.983845478*10^6*x20-5.221405977,
 x26*(-0.3860115660e-1*x14-0.1876793978e-1*x8+0.7836678184e-1)+2.040147478*10^6*x21-10.82168541,
 x26*(-0.1876794098e-1*x14-0.9892449327e-2*x8+0.3810204607e-1)+1.983845478*10^6*x21-8.751240594,
 x15-0.7006679273e-1*x12-.2484955248*x6-1.002451672,
 x16-0.7006679273e-1*x13-.2484955248*x7-1.002451672,
 x17-0.7006679273e-1*x14-.2484955248*x8-1.002451672,
 x18-0.7006679273e-1*x12+.1803623678*x6-1.002451672,
 x19-0.7006679273e-1*x13+.1803623678*x7-1.002451672,
 x2-0.7006679273e-1*x14+.1803623678*x8-1.002451672]:

To solve it with fsolve, Digits must be increased (the system seems to be ill-conditioned).

Digits:=20:
fsolve(sysX,indets(sysX));
{x1 = 1.0112518604279113366, x10 = .50079130545073703755, x11 = -0.42400685353282833597e-1, x12 = 1.0112518604279113366, x13 = .87282458350900007716, x14 = .26867164104671288780, x15 = 1.0733068465024293527, x16 = .96853537340626784353, x17 = .94178275674200828459, x18 = 1.0733068465024293527, x19 = 1.1326128306517202960, x2 = 1.0789746676418324408, x20 = 0.17374637467959646823e-5, x21 = 0.35058655890000000000e-5, x22 = 0.53043642810000000000e-5, x23 = 0.17374637467959646823e-5, x24 = 90.123721949337849443, x25 = 69.574517142133632373, x26 = 49.584072102670758903, x3 = .51046055497717429906, x4 = .91522526886228291076, x5 = .26867164104671288780, x6 = 0., x7 = -.38259166982030833324, x8 = -.31990063204406129237, x9 = 8.4475741104335797627*10^(-7)}

Note that being "almost" linear, using solve instead of fsolve is not a bad idea.

Q:=proc(a,b,c,d,x,y) local t;
   x^2-x+y^2-5*y - evalf(Int(5*exp(sin(t+x)),t=0..b)) end:

Optimization:-NLPSolve('Q(1, Pi/2, 5, 2, x,y)', x= 0 .. 5, y = 0 .. 4);
       [-25.8258678143110600, [x = .750451348084492, y = 2.50000000432642]]

In Maple 2016,

N  := 2:
F  := add(P||k * x^k, k=0..N):
G := [seq(P||k=-1.0 .. 1.0, k=0..N)]:
Explore(plot(F, x=0..1), parameters=G);

works as expected. What version do you use?

Also,

F := a*x:
Explore(plot(F, x=0..1), parameters = [a=0..1]);

works, after replacing (of course) F = a*x with F:= a*x.

But note that you should use [a=0.0 .. 1.0]  otherwise a will take only the two integer values 0 and 1.