[Fulita]

My problem is how should the specific equation be derived using boundary condition. I have the final equation (which is based on Langmuir isotherm) and I know the boundary condition as well so I just need to know how should I reach to the final equation using boundary condition? Step by step.
The model is on the assumptions that (i) in a As–HA binary system, sorbent binding sites form complexes represented as SA when arsenic binding occurs and complexes denoted SH when HA is bound and (ii) equilibrium can be expressed as S+A ↔ SA, where K1 = [SA]/[S][A], S + H ↔ SH, where K2 =[SH]/[S][H] and A+H ↔ AH where K3=[AH]/[A][H]. In this model, considering an interactions resulting in SAH complexes, the equilibrium can be expressed as SA+H ↔ SAH, where K1,2=[SAH]/[SA][H], SH+A ↔ SHA, where K2,1=[SHA]/[SH][A] and S+AH ↔ SAH where K3,3=[SAH]/[S][AH]. Then one can assume S+A+H ↔ SAH and S+H+A ↔ SHA, where K=K1K1,2 = K2K2,1= K3K3,3. Assuming that the sorption system is in equilibrium (there are no net changes of [SA], [SH] and [SAH] with respect to time), the following can be written: d[SA]/dt = 0 , d[SH]/dt=0 and d[SAH]/dt=0 and [S0]=[S]+[SA]+[SH]+[SAH]
In such case, the model can be presented to the form:
Q_As=(qmax Ce [As]{1+(K1/K )Ce [HA] })/({K1+Ce [As]+(K1/K2 ) Ce [HA]+2(K1/K) Ce [As] Ce [HA]})
So I need to know how I should reach to the above final equation using the mentioned assumptions?