The question of self-adjointness is quite often all about the boundary conditions. In order to get the domains of the operator and its adjoint to match, boundary conditions need to be 'distributed' equally between the two.
In this particular case, the first step is to choose a domain $S$ such that the operator $M$ is symmetric, i.e. $(M(u),v)=(u,M(v))$ for all $u,v\in S$.
We have
$$M(u)=\big[\sin^2(x)+\sin(x)-x\cos(x)\big]u-x(xu)'',$$
and when we integrate by parts, we obtain
$$(M(u),v)=(u,M(v))+v'(1)u(1)-v(1)u'(1).$$
These boundary terms have to vanish for $M$ to be a symmetric operator.
Here you already witness the problem of distributing boundary conditions on $u$ and $v$. If we pick $S$ such that $u(1)=u'(1)=0$ for $u\in S$, we only need that $v(1)$ and $v'(1)$ are finite. It seems that choosing $u'(1)=\alpha u(1)$ and $v'(1)=\alpha v(1)$ for some $\alpha\in\mathbb{R}$ is the way to go.
How do we make this rigorous? We typically choose an initial domain $S$ that is a dense subspace of a Hilbert space. In this example, a natural choice for the Hilbert space is $L^2([0,1])$ and a dense subspace on which $M$ can be defined and is symmetric is the subspace
$$S=\big\{u\in W^{2,2}([0,1])\;|\;u(1)=u'(1)=0, \; |u(0)|<\infty \text{, and } |u'(0)|<\infty \big\}$$
of the Sobolev space $W^{2,2}([0,1])$.
With this choice, the domain of $M^{\ast}$ is
$$\mathcal{D}(M^{\ast})=\big\{v\in W^{2,2}([0,1])\; | \; \max\{|u(0)|,|u'(0)|,|v(1)|,|v'(1)|\}<\infty\big\}.$$
The next step is to find the deficiency subspaces $\mathcal{K}_{\pm}\doteq \ker(i\mp M^{\ast})$ and deficiency indices $n_{\pm}=\dim \mathcal{K}_{\pm}$ and then apply theorem $X.2$ (due to von Neumann I believe) on p.140 of Reed, M., Simon, B.: Methods of Modern Mathematical Physics. II. Fourier
Analysis, Self-Adjointness. Academic Press Inc., New York (1975), which I quote here for convenience:
Let $A$ be a closed symmetric operator. The closed symmetric extensions of $A$ are in one-to-one correspondence with the set of partial isometries (in the usual inner product) of $\mathcal{K}_+$ into $\mathcal{K}_-$.* If $U$ is such an isometry with initial space $I(U)\subset\mathcal{K}_+$, then the corresponding closed symmetric extension $A_U$ has domain
$$\mathcal{D}(A_U)=\big\{\phi+\phi_++U\phi_+\big|\phi\in\mathcal{D}(A),\phi_+\in I(U)\big\}$$
and
$$A_U(\phi+\phi_++U\phi_+)=A\phi+i\phi_+-iU\phi_+.$$
If $\dim I(U)<\infty$, the deficiency indices of $A_U$ are
$$n_{\pm}(A_U)=n_{\pm}(A)-\dim I(U).$$
*$\mathcal{K}_{\pm}\doteq \ker(i\mp A^{\ast})$
In your case, most likely you will find $n_+=n_-> 0$ and the above theorem will give you all the possible ways to extend the operator $M$ (i.e. its domain) such that it becomes self-adjoint. I expect (but haven't checked) that you find a one-parameter family of operators corresponding to the boundary condition $u'(1)=\alpha u(1)$.
There are some examples after the quoted theorem in the above-quoted book, and maybe this is already sufficient to get you going. Otherwise, please ask and I will try to give further details.
