Unfortunately, there exist no closed form for this.
The Lagrangean reads
$$
L(w,\lambda)=w^TM_3(w\otimes w)-\lambda(w^T\mathbf{1}-1)
$$
with first order conditions
$$
\begin{align}
\frac{\partial L }{\partial w_i}&=3w^TM_{3,i}w-\lambda \quad \forall i \\
\frac{\partial L }{\partial \lambda}&=w^T\mathbf{1}-1
\end{align}
$$
where $M_{3,i}$ is the $i$th matrix component of the $3$-dimensional skewness tensor. The derivative of $w^TM_3(w\otimes w)$ with respect to $w_i$ is easily verified algebraically, and comparison to a quadratic form.
Effectively, this is a system of quadratic forms:
$$
\begin{align}
w^TM_{3,1}w&=\lambda\\
w^TM_{3,2}w&=\lambda\\
\ldots&=\lambda\\
w^TM_{3,N}w&=\lambda\\
w^T\mathbf{1}&=1
\end{align}
$$
There exist no closed-form solution for this. You could try to solve this equation system using a multivariate Newton Raphson scheme and careful selection of starting values.
Answering your comment:
.... Since there is no closed-form solution for the max skewness portfolio, does that mean that we cannot derive a proof that the max skewness portfolio has higher skewness than the most skewed asset?
At least anecdotically, it is quite easy to show that for a two-asset portfolio, the boundedness of the portfolio skewness is driven by the level of the co-skewness.
Please find below two graphs for a two asset portfolio. In each case, the assets of unit variance, no covariance, and skewness of $S_{111}=0.05$, $S_{222}=-0.05$. In the first graph, the co-skews $S_{112}=S_{122}=0.0$, in the second graph they are $+0.1$ and $-0.1$, respectively. The $x$-axis shows the portfolio weight on asset 1.

As you can see, the question whether or not portfolio skew is bounded by the asset skews is driven by *co-skewness. Again, diversification is the key.
HTH?