The error message is

Error, (in pdsolve/sys/info) the HINTs programmed for PDE systems are: {"as is", `*`, `+`, TWS, tws}; received: HINT = F[1](t)*F[2](r*sin(phi))

i.e., an arbitrary function hint is not available for a system. In the case of a single pde, F(t)*G(r)*H(phi) should work.

HINT=`*` works here to give a separated solution.

If you know it is periodic with period 12, then a Fourier series is the simplest choice. I went to n=3, with 7 parameters for 12 data points, and I guess you could go higher, but at some point you are overfitting.

Download Fit.mw

Your immediate problem is the N_y in the second for statement has an extra space between N and _y. Removing the space gets rid of the error message.

Use KroneckerDelta instead of kroneckerdelta, and since it is in the Physics package, you will have to use with(Physics) at the beginning (or perhaps you don't want it to evaluate?)

You have a matrix element in the quantum mechanics sense, but how do you want to generate a Maple Matrix? Do you want Nx*Ny rows and Nx*Ny columns each with order l=1, p=2; l=1, p=2; l=1, ..., p=N_x, l=2, p=1, etc? So then the  <l,p|H|l_prim,p_prim> element would go in H[l*N_y+p, l_prim*N_y+p_prim] entry (that's probably not right but I'll let you figure it out)

Your if statements don't test l_prim or p_prim.

You can test for even and odd by: if p::even then ...   ; if p::odd then

My tendency would be to just make four nested loops for l, p, l_prim, and p_prim and then just have a single complicated if statement to deal with the various cases.

(With this complex code I'd use 1-D math or a code edit region for ease of debugging.)
 

A degree 6 polynomial doesn't have a general symbolic solution, so you will need to provide values for at least some variables in the coefficients to get an explicit numerical or symbolic solution. (For simpler ploynomials you can give solve the option explicit.)

For future reference, you should upload your worksheet using the big green up arrow.

You are giving dsolve 60 differential equations and also 15 algebraic equations (final[k]), but then you are giving 75 initial conditions. The algebraic equations are just giving values to variables, so you should just substitute these into the differential equations and give dsolve 60 differential equations and 60 initial conditions in 60 variables.

[Edit: OPs uploaded worksheet is modified and is now producing output]

Probably you intended to give kappa and lambda numerical values and forgot. (If you really want them as parameters then use the parameters= option to dsolve - see the help page)

Do(%TextArea0(Enabled)) returns true (and not "true"). I find the Do mechanism easier to integrate into other code for this reason.

Seems there should be a simpler way, but you can do it step by step.

Download sums.mw

sigma[CNT] = 0.1e7; should be sigma[CNT] := 0.1e7;

In general it is helpful for us to solve your problems if you upload your worksheet using the fat green up arrow in the editor.

Since I needed to figure out good guesses and you formulated the ionic strength in molality and the others in moles I got frustrated and cheated by refomulating the whole problem in molalities, which is how I would have set up this problem anyway. These equilibrum problems almost always need high Digits and good guesses. You will see some chemical intuition was also required to solve it.

[Edit: at the end I check it satisfies the OPs eqns]

Download Chem.mw

There is no difference in the singularity if the columns are permuted, so choosing the entries of the first row without considering different orders reduces the number of ways by n! for an n x n matrix.
 

Download HowManyMatricesAreSingular.mw

T is already a Matrix, and stripping off the first two rows and columns can be done with T[3..,3..]

Download CommutatorExample.mw

You need normalization=full, and then, for some reason, another factor of 2 (Every time I do a DFT in a different piece of software, I always need to do a sample to figure this out).

FFT.mw

I couldn't see why yours didn't work, but here's an example that does. The button code is

use DocumentTools in 
nrows:=nrows+1;
A:=Matrix(nrows,2,A);
A[nrows,1]:=0;
A[nrows,2]:=1;
Do(%DataTable0(VisibleRows)=nrows);
end use; 

Download Morerows.mw

If your change ln(2) to ln(2.) then that will force the answer into a single numerical floating point value.