I found the following with the help of Mathematica.
Ansatz:
q_i^* = k_{f,i} c_{p,i} \left( \sum_j a_{i,j} c_{p,j} \right)^{n_i - 1}
p_{i,j}=b_{i,j,1} c_j+b_{i,j,2} c_j^2
Conditions:
p_i(c_j=T)=q_i^*(c_j=T)
\frac{\partial p_i}{\partial c_j}(c_j=T)=\frac{\partial q_i^*}{\partial c_j}(c_j=T)
Coefficients for i=j:
b_{i,i,1}= k_f \left((2-n_i) a_{i,i} T + \sum_{k \neq i}{a_{i,k} c_k}\right) \left(a_{i,i} T + \sum_{k \neq i}{a_{i,k} c_k}\right)^{n_i-2}
b_{i,i,2}=k_f (n_i-1) a_{i,i} \left(a_{i,i} T+ \sum_{k \neq i}{a_{i,k} c_k}\right)^{n_i-2}
Coefficients for i \neq j:
b_{i,j,1}=k_f\frac{c_i}{T} \left((3-n_i) a_{i,j} T+2\sum_{k \neq i}{a_{i,k} c_k}\right) \left(a_{i,j} T+\sum_{k \neq i}{a_{i,k} c_k}\right)^{n_i-2}
b_{i,j,2}= k_f\frac{c_i}{T^2} \left((n_i-2) a_{i,j} T+\sum_{k \neq i}{a_{i,k} c_k}\right) \left(a_{i,j} T+ \sum_{k \neq i}{a_{i,k} c_k}\right)^{n_i-2}
In case c_j<T for more than one index j, I would propose a weighted sum of the corresponding approximations:
p_i^*=\frac{\sum_j{w_j p_{i,j}}}{\sum_j{w_j}}
with weights:
w_j=\max{(0,\frac{T-c_j}{T})}
Please double check for flaws and typos.