The fact is that pdsolve is not very competent with initial conditions
and we must do it ad-hoc:

ds:=pdsolve({sys}) ;
           {x(t, s) = _F1(s) sin(t) + _F2(s) cos(t),
             y(t, s) = _F1(s) cos(t) - _F2(s) sin(t)}
ds0:=eval(ds,t=0);
              {x(0, s) = _F2(s), y(0, s) = _F1(s)}
subs(ds0, [ic]);
                 [_F2(s) = a(s), _F1(s) = b(s)]
eval(ds,%);
             {x(t, s) = b(s) sin(t) + a(s) cos(t),
               y(t, s) = b(s) cos(t) - a(s) sin(t)}

For these kinds of maniputations, the vectors, lists, matrices etc were invented.
Example using column vectors:

gamma__w :=10:
z := <2, 1, 0, -1, -4, -7>:
h := <2, 1, 0,  0,  0,  0>:

gamma__w * (h-z);
                   

You have a polynomial systems with many parameters.
The result of solve would be huge and in terms of RootOfs.

Just to make an idea, let's eliminate C4 between the first two equations:

poly:=(numer@lhs)~([eq1,eq2,eq3,eq4]):
resultant(poly[1], poly[2], C4);

You will see that even this first step gives a huge polynomial.

For numerical results, give values to parameters and use fsolve.

The residuals are much smaller using LSSolve:

eq:=[seq(lhs(EQQ[i]),i=1..36)]:
Digits:=15:
s:=Optimization:-LSSolve(evalf(eq),iterationlimit=100000);

s := [0.960968723053680e-3, [AA[1] = -.347787136052672, AA[2] = 0.510086502851874e-1, AA[3] = 1.16178642934306, AA[4] = .283780970480705, AA[5] = 0.639701227912076e-1, AA[6] = 0.167348852743906e-2, AA[7] = .970111340350232, AA[8] = 0.123075237465221e-1, AA[9] = -2.56229943109674, AA[10] = .119761186424375, AA[11] = 0.514762128059566e-1, AA[12] = 0.653686059489268e-2, AA[13] = 1.17224901537396, AA[14] = 0.123075188025341e-1, AA[15] = -3.10469172766768, AA[16] = .119761196994863, AA[17] = 0.514762189628706e-1, AA[18] = 0.653685044254999e-2, AA[19] = .970010263800571, AA[20] = 0.122243656427024e-1, AA[21] = -2.56235367446998, AA[22] = .119187577004902, AA[23] = 0.512617524055109e-1, AA[24] = 0.636735006533654e-2, AA[25] = 1.17214793882475, AA[26] = 0.122243606766083e-1, AA[27] = -3.10474597103813, AA[28] = .119187587581231, AA[29] = 0.512617585634991e-1, AA[30] = 0.636733990163037e-2, AA[31] = .516890321334600, AA[32] = 0.510086073253494e-1, AA[33] = -1.15838655820066, AA[34] = .283780982288442, AA[35] = 0.639701289588585e-1, AA[36] = 0.167348506238467e-2]]

Digits:=40:
eval( evalf(eq), s[2]);

Unfortunately LSSolve seems to get stuck for Digits>15.

Because we already know what has to be proved, it could be done like this:

Download pde-bug.mw

Edit.  I made a mistake when solving for _F4.
_F3  must be necessarily 0 (due to the presence of y)
So, actually S2 must be

So, indeed, S00 is not general (for example f = t*u+x*y is not contained in S00), but S1 could be after all correct; it would be interesting to check this by hand.
 

subs[eval](F = ( (x,y)->f(x)+g(y) ), F(x,x) );

(Of course [eval] is not mandatory, just to obtain th result evaluated.)

The problem reduces to find/guess  the next term(s) of the sequence:

1, -1/2, 1/12, -1/72, 1/504, -5/18144, 1/27216, -95/19813248

It should be obvious that without extra information this is impossible. Maybe Nostradamus ...

Edit. Example:
f := convert( series(cos(x)+x^2, x, 14), polynom):
L := [seq(coeff(f, x, n), n = 0 .. 12)];

    L := [1, 0, 1/2, 0, 1/24, 0, -1/720, 0, 1/40320, 0, -1/3628800, 0, 1/479001600]

rec := gfun:-listtorec(L, u(n));    # not the expected one
     rec := [{(n^3-26*n^2-449*n+538)*u(n+1)+(-707*n^3-1737*n^2+1283*n+393)*u(n+3), u(0) = 1, u(1) = 0, u(2) = 1/2}, ogf]

ImportMatrix("http://fs3.fex.net/get/245716150875/11071260/data.txt"):
V:=convert(%,Vector);

You have generated only a small number of partitions: 1401400 (=nops(P)),  the total number being bell(N-1) = 27644437.

N is too large for a brute force solution (and with your attempt to insert also the permutations, the needed memory is huge).
I think that a very good (or a heuristic) algorithm is needed here.

Optimization:-Maximize((1-y^2)/(x^2),{x^2 + y^2 <= 1, x>=1/2}, x=1/2..1,  y=-1..1  );
            [ 4., [x = 0.5, y = 0.0] ]

The answer should be obvious without Maple.
For f(u,v) = 1/(u+v)^2+4*u*v-1,

f(0+, 2) = - 3/4 < 0,  f(1, 1) > 0.
f being continuous in the connected domain (0,oo) x (0, oo),  f must be 0 somewhere.

inequal can be used only for two variables.

with(PolyhedralSets):
ps := PolyhedralSet({-5 <= 4*y1+2*y2+3*y3, -2 <= 3*y1+5*y2+2*y3, -1 <= y1+2*y2+y3, 0 <= y2, y1 <= 0}):
Plot(ps);

G:=convert(F,RootOf):
evala(AFactors(G));
convert(%,radical); #optional

#   GroupTheory:-PermParity(Perm(L));
GroupTheory:-PermParity(Perm([1,2,3,5,4]));
                               -1