@Andriy 

I have decided that it is a typo as I can't understand this construction:

N:= (i,j,sigma,f)-> sum(sum(sum(nn(i,j,sigma), sigma= 1..2), j= 1..3), i= 1..f):

You're right: It should be N:= f-> .... What I wrote is nonsense! Sorry about that. Try this:

nn:= ()-> `if`([args]::list(integer), ap.am(args), 'procname'(args)):

N:= proc(f)
local i,j,sigma;
     Sum(Sum(Sum(nn(i,j,sigma), sigma= 1..2), j= 1..3), i= 1..f)
end proc:

I can't help with the simplification of the nn operators.

@Andriy 

I have decided that it is a typo as I can't understand this construction:

N:= (i,j,sigma,f)-> sum(sum(sum(nn(i,j,sigma), sigma= 1..2), j= 1..3), i= 1..f):

You're right: It should be N:= f-> .... What I wrote is nonsense! Sorry about that. Try this:

nn:= ()-> `if`([args]::list(integer), ap.am(args), 'procname'(args)):

N:= proc(f)
local i,j,sigma;
     Sum(Sum(Sum(nn(i,j,sigma), sigma= 1..2), j= 1..3), i= 1..f)
end proc:

I can't help with the simplification of the nn operators.

@Markiyan Hirnyk asked: Could you explain the occurrence of or in your code?

I also wondered at it at first. But pay close attention to the direction of the inequalities that are connected by or. Preben had to do it this way because two-argument arctan always returns a value between -Pi and Pi. So the shaded region is considered as two regions by inequal.

@Markiyan Hirnyk asked: Could you explain the occurrence of or in your code?

I also wondered at it at first. But pay close attention to the direction of the inequalities that are connected by or. Preben had to do it this way because two-argument arctan always returns a value between -Pi and Pi. So the shaded region is considered as two regions by inequal.

@Jimmy I think (just guessing here) that when the fit is very good, as it is in your situation, that it can work without the parameter ranges.

@Jimmy I think (just guessing here) that when the fit is very good, as it is in your situation, that it can work without the parameter ranges.

@Andriy But you never tested what I actually wrote! You have N:= f-> ..., but I wrote N:= (i,j,sigma,f)-> ....

@Andriy But you never tested what I actually wrote! You have N:= f-> ..., but I wrote N:= (i,j,sigma,f)-> ....

Your procedure can be much simpler:

P:= proc(N1,N2)
local q,r;
     if N1 > N2 then error end if;
     q:= iquo(floor(N2)-ceil(N1)+1, 12, 'r');
     2*q+`if`(r>3, `if`(r>8, 2, 1), 0)
end proc;

Your procedure can be much simpler:

P:= proc(N1,N2)
local q,r;
     if N1 > N2 then error end if;
     q:= iquo(floor(N2)-ceil(N1)+1, 12, 'r');
     2*q+`if`(r>3, `if`(r>8, 2, 1), 0)
end proc;

@mehdi jafari I think that one problem is that in your eq1 (in problem 3.mws) there occurs both a[1](t) and a[1] without the (t). Is that what you really intended? I can't think of any reason why someone would want to use the same variable as both a function and a non-function in one expression.

For your sample probem with eq1, eq2, and eq3, can you write down "by hand" what the Matrix A should be? That is, does there exist A such that A . < a1(t),a2(t),a3(t) > = < eq1,eq2,eq3 >. It is not clear to me that such A exists.

@mehdi jafari I think that one problem is that in your eq1 (in problem 3.mws) there occurs both a[1](t) and a[1] without the (t). Is that what you really intended? I can't think of any reason why someone would want to use the same variable as both a function and a non-function in one expression.

For your sample probem with eq1, eq2, and eq3, can you write down "by hand" what the Matrix A should be? That is, does there exist A such that A . < a1(t),a2(t),a3(t) > = < eq1,eq2,eq3 >. It is not clear to me that such A exists.

Another way is to wrap the sum in a procedure. I think that that is more robust.

I see that you also asked for suggestions on simplifying your code. Well, all your procedures can be made one-liners with arrow -> notation. nn can be simplified further because function application distributes over operators. I also changed the grammar of your comments from passive declarative to active imperative.

restart:
with(Physics):
Setup(mathematicalnotation= true):
Setup(anticommutativeprefix= Psi):

#Define the annihilation operators, which depend on indicies (i, j, sigma).
am:= (i,j,sigma)-> Annihilation(Psi, 6*(i-1)+2*(j-1)+sigma, notation= explicit):

#Define the creation operators, which depend on indicies (i, j, sigma).
ap:= (i,j,sigma)-> Creation(Psi, 6*(i-1)+2*(j-1)+sigma, notation= explicit):

#Program the occupation number at state (i, j, sigma).
nn:= ap.am:

#Define an operator for the total number of particles.
N:= (i,j,Sigma,f)-> sum(sum(sum(nn(i,j,sigma), sigma= 1..2), j= 1..3), i= 1..f):

Another way is to wrap the sum in a procedure. I think that that is more robust.

I see that you also asked for suggestions on simplifying your code. Well, all your procedures can be made one-liners with arrow -> notation. nn can be simplified further because function application distributes over operators. I also changed the grammar of your comments from passive declarative to active imperative.

restart:
with(Physics):
Setup(mathematicalnotation= true):
Setup(anticommutativeprefix= Psi):

#Define the annihilation operators, which depend on indicies (i, j, sigma).
am:= (i,j,sigma)-> Annihilation(Psi, 6*(i-1)+2*(j-1)+sigma, notation= explicit):

#Define the creation operators, which depend on indicies (i, j, sigma).
ap:= (i,j,sigma)-> Creation(Psi, 6*(i-1)+2*(j-1)+sigma, notation= explicit):

#Program the occupation number at state (i, j, sigma).
nn:= ap.am:

#Define an operator for the total number of particles.
N:= (i,j,Sigma,f)-> sum(sum(sum(nn(i,j,sigma), sigma= 1..2), j= 1..3), i= 1..f):

I answered your followup question in the separate thread.

Don't crosspost, i.e., don't post the same thing in multiple places. I deleted three (!!!) copies of your followup question from this thread.