I don't think that assuming works with fsolve, but evalc will convert the equations into something that fsolve will handle. In particular, it replaces abs with an equivalent sqrt expression and replaces argument with an equivalent arctan expression.

fsolve(evalc~({eq1, eq2}), {a, z});

     {a = 0.561266211604431e-1, z = .313474043774231}

If x is a decimal number, then the number of places after the decimal point is -op(2,x). This is not necessarily the same as the number of digits of precision of x.

You shouldn't use the same symbol, x, to represent both the unknown in your expression and the initial value of that unknown. I prefer to use x0 to represent the initial value. Here's some corrected code:

simple_iterations:= proc(
   fun::algebraic, x0::complexcons, 
   {epsilon::positive:= 1e-5}, {itermax::posint:= 20000}
)
local N:= indets(fun, And(name, Not(constant))), f, iter, rez, x;
   if nops(N)>1 then error "Too many variables." end if;
   f:= evalf@unapply(fun, [N[]]);
   rez:= x0;
   x:= f(x0);
   for iter to itermax while abs(rez-x) > epsilon do
      rez:= x;
      x:= f(x)
   end do;
   if iter = itermax+1 then error "Did not converge" end if;
   x
end proc:

simple_iterations(1/sin(Pi*x/180), 7);

I only have time for a quick Answer right now. The answer to your second question is a definite "no"; it's typical to just use a "large" value to approximate infinity. But the value that you've used is much too large. The values used are typically 5 or 10.

I'll look at your problem in more detail later. One thing that stands out is that none of your IBCs fix the value of x.

Maple's numerical PDE solver is very limited. There must one spatial (space-based) and one temporal (time-based) independent variable, at least implicitly. Which do you consider the spatial, and which do you consider the temporal?

Here is a highly generalizable process for finding the area enclosed by curves. But you must still look at the plot to decide which boundary is which.
 

Download Area_between.mw

A boundary that is a vertical line becomes one of the limits of integration.

If there is any hidden code associated with worksheet, then it can be accessed from the menu Edit -> Startup Code. Many (maybe most) worksheets have their code in plain view right in the worksheet.

Here is how to classify binary matrices by their rank over GF(2).
 

Download rank_mod_2.mw

If someone knows how to make this more efficient by using ArrayTools:-Alias or by other inplace operations, I would appreciate seeing it.

How about

(T-> (p-> [T[..-2][], p])~(combinat:-permute(T[-1])[..max(2,nops(T[-1])!)/2])[])~(Tours);

But, there must be a better way to your overall goal. If you describe that goal better, maybe I can figure out a better way.

min({L[ ]} minus {0});

I think that you're having a problem that the other respondents aren't because you're saving the module to a library and they aren't. Global variables aren't part of the module and thus won't be saved unless they are explicitly mentioned. There is no need for global variables. Try this:

MonPackageFonctions:= module()
option package;
local GeometricData, r, l, ModuleLoad:= ()-> SetGeometricData(0.2, 1);
export 
   GammaNum:= xx->
      evalf(eval(Pi-arccos((-xx^2+l^2-r^2)/xx/r/2), GeometricData)),
   psiNum:= xx->
      evalf(
         eval(
            -arcsin(sqrt(-(xx^4-2*xx^2*l^2-2*xx^2*r^2+l^4-2*l^2*r^2+r^4)/xx^2/r^2)/2/xx/l),
            GeometricData
         )
      ),
   SetGeometricData:= proc(r::algebraic, l::algebraic)
      GeometricData:= [thismodule:-r= r, thismodule:-l= l];
      NULL
   end proc
;
   ModuleLoad()
end module:

This exports a procedure that lets you set the GeometricData externally but also uses r= 0.2 and l=1 as defaults.

If the file is named myfile.bat then do

system("myfile.bat");

There are many possibilities. The criteria are mostly aesthetic. How about this?
 

Download Tours.mw

In differential equations in Maple, the dependent variables must always be used in functional form, i.e., x(t) and y(t) rather than x and y. Therefore, your code should be changed to

with(inttrans);
eq5 := diff(x(t), t)+diff(y(t), t$2) = 5*exp(2*t);
eq5s := laplace(%, t, s);
eq6 := diff(x(t), t)-x(t)-(diff(y(t), t))+y(t) = 8*exp(2*t);
eq6s := laplace(%, t, s);
solve({eq5s, eq6s}, {laplace(x(t), t, s), laplace(y(t), t, s)});
subs({x(0) = 2, y(0) = 1, (D(y))(0) = 1}, %);
eq3 := invlaplace(%, s, t);

If you make this simple change, then the last line will give you the answer that you're looking for.

A more-sophisticated approach---which might be useful if you want to re-use this operation a lot---is to use the RandomTools package:

RandomTools:-AddFlavor(
   A = module()
   local r:= rand(-90..90);
   export Main:= ()-> (a-> [a, convert(a*degrees, radians)])(r());
   end module
):
randomize():
n:= 4:
RandomTools:-Generate(list(A, n));