kill(all)$

ode : part6,intro: ricatti 
peter.vlasschaert@gmail.com,05/02/2017

# Theory Ricatti : ODE (reaction kinetics)

(dy/dx) = p(x)*y^2+q(x)*y+r(x)
Why: Consecutive Reactions and Ricatti
                      k1   k2
reaction kinetics : A -> B -> C
dCA/dt = - k1*CA^n          (1)
dCB/dt = k1*CA^n-k2*CB^m    (2)
dCC/dt = k2*CC^l            (3)
                                              
integer values :  n,m,l 
concentration  : CA,CB,CC

load('contrib_ode)$
assume(t>0,CA>0,CB>0)$

equation (2) n=2,m=1, by 't=0 -> CA=CA0'

eq1: CA = CA0/(1+k1*CA0*t);
eq2:'diff(CB,t)=k1*CA^2-k2*CB^1;
eq3:lhs(eq2)-part(eq2,2,2)=rhs(eq2)-part(eq2,2,2);

eq4:subst(rhs(eq1),CA,eq3);
eq5: lhs(eq3)=rhs(eq4);

eq6:contrib_ode(eq5,CB,t)[1];
method;

function :build in maxima = 'expintegral_e(r,z)'
info :

https://en.wikipedia.org/wiki/Exponential_integral
http://maxima-online.org/help/index/expintegral_e
http://functions.wolfram.com/GammaBetaErf/ExpIntegralE/
http://encyclopedia2.thefreedictionary.com/Exponential+integral

eq7:expintegral_e(2,-(k2*(CA0*k1*t+1))/(CA0*k1));

find %c,by t=0 and CB=0 in eq6

eq8:ev(eq6,t=0,CB=0);
eq9:solve(eq8,%c)[1];

solution : n=1,m=2 for eq5

eq10:ev(eq6,eq9);