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.