I would like some tips in figuring out a good algorithm to find the solution of the following system. Let $\theta$ be a constant in $(0,1)$, let $i,l=1,...,N$, let $a_{l}$ and $b_{i,l}$ be some constants with $a_l \in (0,1)$ and $b_{i,l}>0$ for all $i,l$. Consider the following system of equations in unknowns $z_{1},...,z_{N}$ and $p_{1},...,p_{N}$:$$a_{l}=\sum_{i}\left(\frac{b_{i,l}}{\sum_{i}b_{i,l}}\right)\frac{1}{1+p_{l}z_{i}},l=1,...,N$$ $$\frac{1}{1+z_{i}^{\theta}}=\sum_{l}\left(\frac{b_{i,l}}{\sum_{l}b_{i,l}}\right)\frac{1}{1+p_{l}z_{i}},i=1,...,N$$
I am focusing on real and positive solutions.
Newton methods require lots of guidance to get good initial values, so I have focused on iteration algorithms. One can construct a mapping $ p \rightarrow p'$ as follows: first take $p$ and solve for $z$ from the second set of equations and then go to the first set of equations with that $z$ and solve for $p'$. If $\theta<0$ then this is a contraction mapping, and it stays in some box, implying that there is a unique solution and one can find it using this iteration. But now I want to focus on $\theta\in(0,1)$.
One idea is to again construct an iteration from $p$ to $p'$ but not reversing the one above: first take $p$ and solve for $z$ from the first set of equations and then go to the second set of equations with that $z$ and solve for $p'$. But here solving each set of equations is more involved as I cannot solve them separately, and so again I am at the mercy of fsolve which requires extra work to get appropriate initial values, making the algorithm too slow. It may also be that for some $p$ the first equation has no solution.
Perhaps better to use the second equation to solve for $z_i$ for a given $p$ and the first equation to solve for $p_l$ for a given $z$. The question is how to combine these two functions $z_i(p),p_l(z)$ into a robust and quick algorithm to get to the solution.
Any tips would be much appreciated.