I simplified my setup as much as possible. Please check lambdas.mw.
While I think I managed to obtain some analytical solutions, they look a bit strange for two reasons:
1) They do not depend on the exogenous parameters as I expected. In fact, mu_jk and mu_ki should only depend on q_0jk and q_0ki, while lambda_jk and lambda_ki should only depend on BigSigma_0jk, BigSigma_0ki, smallsigma_ujk and smallsigma_uki.
2) Strong dependence on q_0jk and q_0ki: if I were to setup these two parameters to zero or to the same value I can't obtain solutions anymore (especially for the lambdas). Does it mean that they are not really "free" parameters?
I noticed that if I combine the two equations from the FOCs of mu_jk and mu_ki into one system (is this even legit?), I get q_0jk = - q_0ki * (lambda_jk / lambda_ki). This is also easy to see if I apply the calibration at the beginning of the script (remove hashtags on all the params with the exception of q_0jk and q_0ki) and then divide lambda_jk by lambda_ki. Why?
I am quite sure that the computations are correct (I checked multiple times), but I am now questioning my setup. In which ways does my setup differ from the one below?
Essentially, I am trying to extend the following problem. As you see below, mu depends only on p_0 (the one-dimensional equivalent of my q_0jk and q_0ki) and lambda depends only on BigSigma_0 and smallsigma_u (the one-dimensional equivalents of my BigSigma_0jk, BigSigma_0ki, smallsigma_ujk and smallsigma_uki).