@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.

@Jimmy Here is your file with my modifications. It is essentially the same code as the one that I did with fake data. I needed to expand your parameter ranges in order to get a good fit. I hope the new parameter values make sense for your problem.

Download TE.mw

I don't think it fair that you criticize my answer so harshly when you didn't post your code. If I had seen that you were using big integers, I would've figured out the problem.