It's easy to calculate the modulus.  It is also easy to verify that it is not always less than 1.  Why do you expect it to be less than 1?

z := coth(x+I*y);

Compute the (square of the) modulus:

z*conjugate(z):
expand(%) assuming x::real, y::real:
m := simplify(%);

Now m is the square of the modulus of z at any x and y.  Let's evaluate it at a point on the boundary of your square domain with n=0:

subs(x=(n+1/2)*Pi, y=0, n=0, m):
evalf(%);

The result is greater than 1.

Try with n=1:

subs(x=(n+1/2)*Pi, y=0, n=1, m):
evalf(%);

Again, the result is greater than 1.

At any given location the daylight fraction varies with the time of the year.  Days are longer in the summer and shorter in the winter.  The calculations in your worksheet do not account for that.

I have attached a worksheet which introduces a function Z(tau, lambda, t), where tau is the tilt of Earth's axis, lambda is the latitude of a location, and t is the time of the year.  In a full year t goes from 0 to 2*Pi.  The function Z is the fraction of daylight at a latitude lambda at time t of the year.

Here are the graphs of Z(Pi/6, Pi/4, t) and Z(Pi/6, lambda, Pi/2):

daylight2.mw

Note added on Dec 28:

I wasn't happy with the ugly calculations in the worksheet which I posted above.  I examined the details and saw that the complications were due to a non-optimal choice of basis vectors.  With a better choice of basis vectors we get a cleaner calculation.  See:

daylight3.mw

The differences between the two worksheets are cosmetic.  The resulting graphs are identical.

Perhaps this is what you are looking for:

http://www.mapleprimes.com/questions/208795-Periodic-Piecewise-Function

You may refer to the "normal" sum through :-sum, as in:

f := t,n -> :-sum(MyFunction(k,t),k=1..n);

Independently of this issue, you should parenthesize the "t,n" part:

f := (t,n) -> :-sum(MyFunction(k,t),k=1..n);

Optimization:-LSSolve gives a reasonably good answer:

PA = 0.279246445735439
PB = 0.072848034900520
PO = 0.648114970425358

See the worksheet for details.

mw.mw

There are very many ways of doing this.   Here is one:

    eval(D(y)(x)=sin(x-y(x)), y=(x->a*x+b));

Beware that your instructor may not be amused by a solution that goes beyond what was taught in class.

In ode1 you have a factor:

At z=0 you are setting T(0)=350, whereby the term above evaluates to

which is complex-valued.  Do you really want complex-valued solutions?

Code Edit Region, which you will find it under the Insert menu, may be of interest to you.  It inserts a editable block into the worksheet, within which the code that you type is autoindented.  (Well, forward indentation is automatic; for backward indentation hit the Backspace key.)  There are quite a few additional features as well; see the help page on Code Edit Region.  Hit Ctrl-E to execute the block's code.

plots:-implicitplot3d([y=x^2, x^2+4*y^2+4*z^2=16],
    x=-4..4, y=-3..3, z=-3..3, scaling=constrained,
    grid=[50,50,50], style=surface);

Alternatively, with parametrized surfaces:

plots:-display([
    plot3d([4*sin(t)*cos(s), 2*sin(t)*sin(s), 2*cos(t)], t=0..Pi, s=0..2*Pi),
    plot3d([x,x^2,z], x=-2..2, z=-3..3)
], scaling=constrained);

restart;
u := Heaviside;
f := t -> u(t);
g := t -> 2*((t-1)*u(t-1) - (t-4)*u(t-4));
int(f(t-q)*g(q), q=-2..10) assuming t>-2, t<10;

PS: You posted a picture of your code, therefore I had to retype the code myself which is a bother.  In consideration for people who may want to help you, post real code that can be copy/pasted, like what I have done above.

params := p=1, w[u]=0.8, w=0.5, v=0.5, D[A]=1000, D[B]=1000, H=2500;
plot3d(eval([Profit_A, Profit_B], [params]), P[A]=0..P[B]-0.01, P[B]=0..0.5);

Your eq₅ is:

eq₅ := diff(DD(t),t) = d₁ E(t) + d₂ In(t);

Introduce a new unknown, let's call it C(t), so that the derivative of C(t) equals DD(t), and therefore C(t) is the antiderivative of DD(t).  This appends the equation:

eq₆ := diff(C(t),t) = DD(t);

to your system.  For the initial condition take Init₆ := C(0) = 0. Then solve the system of 6 equations in 6 unknowns, and plot C(t) as you wish.

I assume that you are loading the PDEtools package because you intend to solve a PDE.  Let mypde be that PDE.  Then instead of solving mypde,  solve the system consisting of { mypde, diff(u,x)=0 }.  This will force diff(u,x) and all its higher order derivatives to be zero throughout.

The short answer to your question is:

q1 := eval(q, omega=2646.408147):
NullSpace(q1);

That said, I should add that I wouldn't put trust in this result.  I fact, you shouldn't put trust in your own matrix q, which contains numerical entries like 10^(-118)  and  10^17.  To put things in perspective, let's note that the distance to the Moon is "only" 3.844 x 10^5 km.  Note the 10^5.  Distance to the Sun is "only" 1.52 x 10^8 km.  Distance to the nearest star is 4 x 10^13 km.  What does your 10^17 measure?  And in what units?

I suggest that instead of continuing with the current calculation (which I consider insane) go back to the derivation of the equations and see where those numbers come from.  If you were doing astronomy, for example, you wouldn't express distance to Neptune in millimeters.  Instead, you would say that the distance to Nepute is 30 AU.  One AU (Astronomical Unit) is the distance between the Earth and the Sun.

Huge variations in numerical coefficients result in numerical instability and nonsensical results.  Choose your units correctly to obtain reliable answers.

Let's see what ZZ looks like as a function of omega:

plot(ZZ, omega=0..10);

We see that the graph confirms your calculation that the positive roots are at:

     1.360526281, 4.077926876, 6.798915756, 9.453815360, 9.864454570

Now, you plug the root 1.360526281 in ZZ and you get something like -1.11*10^58 rather than zero.  There is nothing wrong in that.  Considering that the function's amplitude is of the order of 10^67, a number like 10^58 is one billionth of the amplitude, and therefore it is zero for all practical purposes.