I was given that the free energy per lattice site of the 3 dimensional-dimensional 3 state-state Potts model in the mean field approximation is $$f\equiv \frac{F}{\text{# sites}} = \sum_{k=0}^{2}(-3Kx_k^2 + \frac{1}{\beta}x_k \ln x_k\Big)$$ where here $K$ is from $$H = -K \sum_{\langle ij\rangle} \delta_{\sigma_i\sigma_j}$$ and $x_k$ are the order parameters, the proportion of lattice sites in the state k. To find the temperature where a phase transition occurs, I figured that we take derivatives: $$\frac{\partial }{\partial x_{k_0}}f = -4Kx_{k_0} + \frac{1}{\beta}\ln x_{k_0} + \frac{1}{\beta} = 0\Rightarrow 4K\beta x_{k_0}-\ln x_{k_0}=1$$ but this seems difficult so solve. It's also not clear to me what to do since we have three order parameters.
Context: Consider a system on lattice sites $i \in \mathbb{Z}^d$, where the state $\sigma_i$ of a site can take values in $\{1,\dots, q\}$. We have a Hamiltonian on this system given by $$H = -K \sum_{\langle ij\rangle} \delta_{\sigma_i\sigma_j}$$ where the sum is over pairs $i,j$ of neighboring lattice sites and $K>0$.