k = 5000 Cao = 0.01 alpha = 0.01 epsilon = yao*1 d(p) / d(W) = -alpha*(1+epsilon*X)/2/p p(0) = 1 yao = 0.33 W(0) = 0 Cc = Cao*3*X*p/(1+epsilon*X) thetaB = (1-yao)/yao Cb = Cao*(thetaB-X)*p/(1+epsilon*X) W(f) =80 Ca = Cao*(1-X)*p/(1+epsilon*X) ra = -k*(Ca*Cb -Cc^3/Kc) d(X) / d(W) = -ra/Fa0 X(0) = 0 Fa0 = Cao*vo Kc = 0.05 vo=1500