@sand15: No, phi is well-defined, given by

simplify(convert(eval(phi,k = I),exp));

this being, by the way, the answer to the last question of the original poster. But I guess I know where your mistake/confusion may originate: Note that the expression

which is indeed ill-defined for k = I is not phi itself, but the square of the absolute value of phi evaluated for real-valued k (so k = I cannot even be taken for that expression).

@Zeineb: Assuming that k is real-valued, the square of the absolute value of phi is given by

simplify(evalc(abs(phi))^2);

This should vanish outside some region, then, or at least fall of, say, exponentially fast, for some value(s) of k. I fail to see how that is possible when the graph for tanh(x)^2 is as follows:

Can k be complex-valued? If that is allowed, then letting k = a + b*I, with a and b real-valued, the square of the absolute value of phi becomes instead

k := a + I*b:
normSquare := simplify(evalc(abs(phi)^2));

That is much more promising, for with (a,b) = (0,1), this norm square equals

convert(eval(normSquare,{a = 0,b = 1}),exp);

the graph of which is

So for k = I, phi may indeed be interpreted as a bound state. Whether this is the complete solution for k remains to be investigated.

@Zeineb: Much better. With

Eq := diff(phi(x,k),x$2)+(k^2+2*sech(x)^2)*phi(x,k),

phi as given is indeed a solution. So let us return to your question concerning being a bound state or not: Does bound mean a solution localized to a particular region, possibly with exponential tails dying off ouside such a region, or what?

I am unable to verify that phi as given is a solution to the differential equation:

Eq := diff(phi(x,k),x$2)+(k^2+2*sech(x))*phi(x,k):
simplify(eval(Eq,phi(x,k) = (I*k-tanh(x))*exp(I*k*x)/(I*k-1)));

The expression in paranthesis does not vanish.

As suggested in a previous post of mine (an answer to a previous question of yours on which I would have appreciated some response from you): In order to avoid loosing leading zeros, use lists throughout.

@ecterrab: Thanks for a comprehensive answer. I think I will frequently return to this post of yours, and from its listed four main entry points try to further pave some way forward for me in the Physics environment. Looking forward to an e-textbook, nonetheless :-).

@vv: The underlying O(8) symmetry of the problem will, if not explicitly broken 'by hand', result in any solution being parametrized by 28 continuous parameters, as well as one discrete parameter. That, I guess, is not going to make the problem any easier to either obtain or interpret.

@vv: With your edit, I now comprehend what you have in mind. Do you have any concrete ideas about how to proceed? What explicit construction method, if any, do you have in mind?

@vv: I see what you mean. So my coding really was much ado about nothing, somewhat embarrasingly, I guess.

What precisely do you then mean by "It would be nice to use Maple to find such a configuration"?

@vv: Perhaps I am being stupid here, but I cannot see how you can conclude ||x - y|| >= sqrt(2) from ||x|| = ||y|| = sqrt(2). At the very least, the triangle inequality gives ||x - y|| <= ||x|| + ||y|| = 2*sqrt(2), the inequality sign pointing 'the wrong way'. But perhaps you are (tacitly) using some properties of root/weight space that I am not aware of?

@mskalsi: I am afraid that I will be of little help on these matters, as I have no working experience, mathematically, with quotient algebras. But I have one comment: I think that the algebra you are taking the quotient with respect to have to be an ideal in the overall algebra, i.e., W1 has to be an ideal in SubalgebraNormalizer(W1). Using Query(..."Ideal"), see the help page on Query, I found yesterday that that was not the case. But I, too, may be using things the wrong way.

Have you tried the command QuotientAlgebra, see its help page?

@Muhammad Usman: Oh, I see, I was to quick (did not download your attached worksheet). For a (probable) solution, I happily refer you to the above post by tomleslie.

@mskalsi: You are welcome.

Concerning literature other than text books, I have benefitted from the review article 'Group theory for unified model building', by R. Slansky.

@mskalsi: Concerning your post 'Decomposition':

What you do is essentially the same as I do in my post 'Decomposition' above, but without the last crucial part of finding the new transformed set of generators, the ones in which the Lie algebra decomposes into indecomposable subalgebras.

This last step can be done by multiplying the generators (arranged as a Vector) by the Matrix contained in the quantity 'decomposition', Eq. (5) of yours, i.e., multiplying by decomposition[1]. This is what I do in my post 'Decomposition'. After that, you can initialize a new Lie algebra in terms of these new generators, and then from that print the multiplication table. This latter will differ from your Eq. (4), most importantly by having decomposed the 4x4-block further. But now we know, having used Decompose(), that the algebra is really decomposed into indecomposable subalgebras.

PS: I have trouble copy-pasting multiplication tables from Maple to this text editor so that they render nicely. That is the reason why there were no such in my post 'Decompose' above. How do you produce your images in your post above?