This is a job for fsolve, not solve. You should've uploaded the entire system so that I could test it. You may also try DirectSearch, which needs to be downloaded from the Maple Applications Center.

If you had only allowed non-strict inequalities, then this would simply be a matter of calling Optimization:-Minimize with a trivial objective function. Allowing strict inequalities, it's a little more complicated. Here's a procedure for it:

FeasiblePoint:= proc(
     Cons:= set({`<`, `<=`, `=`}),
     {epsilon::positive:= 10.^(2-Digits), epsilonmax::positive:= 1}
)
local
     newCons:= subsindets(Cons, `<`, ineq-> lhs(ineq) <= rhs(ineq) - epsilon),
     Sol:= Optimization:-Minimize(0, newCons)[2]
;
     if not hastype(Cons, float) then Sol:= convert(Sol, rational) end if;
     if andmap(is, eval(Cons, Sol)) then Sol
     elif epsilon < epsilonmax then thisproc(Cons, 'epsilon'= 10*epsilon)
     else error "No feasible point found; epsilonmax exceeded."
     end if
end proc:   

Example of use:

FeasiblePoint({ x < 1 - y, y < 0, y + 1 < x });

If kx=0 or ky=0 or both, then the integral is 0, as can be determined by direct evaluation. The result returned for the general case is still valid for kx=0 and/or ky=0. So there's no need for piecewise.

The absolute value below is in case you don't know which function is on top. That's trivial in this case, but I wanted to be a bit more general.

F:= (-x^2+9) - (x+3):
abs(int(F, x= `..`(solve(F, x))));

TangentLine:= proc(f::algebraic, at::name= algebraic)
local x, a;
     (x,a):= op(at);
     eval(diff(f,x), x= a)*(x-a) + eval(f,x= a)
end proc:

NormalLine:= proc(f::algebraic, at::name= algebraic)
local x, a;
     (x,a):= op(at);
     eval(f, x= a) - (x-a)/eval(diff(f,x), x= a)
end proc:

f:= ln(3*x)+3; a:= exp(1);
                        f := ln(3 x) + 3
                          a := exp(1)

TL:= TangentLine(f, x= a);

NL:= NormalLine(f, x= a);

plot([f, TL, NL], x= 1/2..10, y= 1/2..10, legend= [f, tangent, normal]);

The following computation strongly suggests that the value that you seek is undefined:

 MG:= MeijerG([[-.3+eps], []], [[.8, 1.3, -.8, -.3, -1.3], []], 1.):
eval(MG, eps= 10^k) $ k= -6..-2;

eval(MG, eps= -10^k) $ k= -6..-2;

Like this:

Ys:= [1,2]: #List all y values here.
Zs:= [3,4]: #List all z values here.

G:= (y,z)-> sqrt(D[1](x)(y,z)^2*dy^2 + D[2](x)(y,z)^2*dz^2):
zip(G, Ys, Zs);

See the command ?map. Example:

V:= <1,2>;

map(tan, V);

What you are trying to do only makes sense to me if E is sorted. So, assuming that E is sorted, a binary search will be much faster in the long run than the linear searching already mentioned. What you want is

1+ListTools:-BinaryPlace(E,x);

Example:

E:= [0,2,7,15,26,40]:
1+ListTools:-BinaryPlace(E,5);
               3

I'll show you a way to get the reduced row echelon form and leave it up to you to interpret it.

A:= < 1,2,3,1; 4,5,6,1; 7,8,9,1 >:

(How did I get that matrix from your equations?)

MTM:-rref(A);

The above command is equivalent to LinearAlgebra:-ReducedRowEchelonForm but is a lot less to type.

Include

with(Optimization):

at the top of your worksheet.

Enclose the _X in single back quote characters, like this: `_X`. On American keyboards, this character is in the upper left corner, under Escape. Anything inside the quotes is treated as a variable, regardless of any other meaning(s) that the characters have.

Use square bracket indexing for subscripting and exponents (via ^) for superscripting. Here's an example:

plot(x^2, x= -1..1, labels= [x[old], x[old]^2], labelfont= [TIMES,ITALIC,16]);

Maple will do the integral easily if you assume x > 0. I also changed the floating-point constants to exact rationals, although this is not strictly necessary.

Int((x-tau)^(1/2-1)*(tau^2+tau^(3/2)*8/3-tau-2*tau^(1/2)), tau= 0..x);

value(%) assuming x > 0;

What do you mean by (x)(x)? Do you mean x times x? That should be either x*x or x^2. Juxtaposition does not imply multiplication. Extraneous parentheses do imply function calls, i.e., (x)(x) is treated as if x were a function with itself as the argument.

The numeric constant 5.6*10^-4 needs to be 5.6*10^(-4) or, better yet, 5.6e-4.

solve(5.6e-4 = x^2/(0.2-x), x);