Laplce transform solution to this system of Nonlinear ODEs I want to solve this system of advective-diffusive-reactive equations analytically:
$$\left(\alpha - k_0c_B\right)c_A+v\frac{dc_A}{dx}-D\frac{d^2c_A}{dx^2} = f_A $$
$$\left(\alpha - k_0c_A\right)c_B+v\frac{dc_B}{dx}-D\frac{d^2c_B}{dx^2} = f_B $$
$$k_0c_Ac_B+\alpha c_C+v\frac{dc_C}{dx}-D\frac{d^2c_C}{dx^2} = f_C $$
Where $\alpha$, $k_0$, $v$, $D$, and $f_{A/B/C}$ are nonzero constants. Here are the boundary conditions:
$$c_A(x=0)=0,\;c_A(x=1)=0$$
$$c_B(x=0)=1,\;c_B(x=1)=0$$
$$c_C(x=0)=0,\;c_C(x=1)=0$$
if $k_0$ were zero, this becomes a system of independent and linear ODEs which I could solve using for example the Laplace Transform. I heard that for nonlinear ODEs, one cannot necessarily use Laplace Transforms without some modification to the original system (like linearization).
If there is a better way to solve the above system, I would like to know how, thanks!
 A: One thing you can try is expand each of $c_A$, $c_B$, $c_C$ in a series in powers of $k_0$.  The equations for each power of $k_0$ will be linear in the coefficients of that power in $c_A$, $c_B$, $c_C$ given the lower-order coefficients.  
EDIT:  To reduce the number of subscripts, I'll write $k_0 = \epsilon$, $c_A = A$, $c_B = B$, $c_C = C$.
Thus  $$\eqalign{A &= A_0 + A_1 \epsilon + A_2 \epsilon^2 + \ldots\cr
B &= B_0 + B_1 \epsilon + B_2 \epsilon^2 + \ldots\cr
C &= C_0 + C_1 \epsilon + C_2 \epsilon^2 + \ldots\cr
A_0(0) &= A_1(0) = A_2(0) = 0,\ A_0(1) = A_1(1) = A_2(1) = 0\cr
B_0(0) &=1, B_1(0) = B_2(0) = 0,\ B_0(1) = B_1(1) = B_2(1) = 0\cr
C_0(0) &= C_1(0) = C_2(0) = 0,\ C_0(1) = 1,\ C_1(1) = C_2(1) = 0\cr
}$$
and the differential equations become
$$ \eqalign{
\alpha\,A_{{0}} \left( x \right) &+v{\frac {\rm d}{{\rm d}x}}A_{{0}}
 \left( x \right) -  D  {\frac {{\rm d}^{2}}{{\rm d}{x}^{
2}}}A_{{0}} \left( x \right) -f_{{A}} = 0\cr
\alpha\,B_{{0}} \left( x \right) &+v{\frac {\rm d}{{\rm d}x}}B_{{0}}
 \left( x \right) -  D {\frac {{\rm d}^{2}}{{\rm d}{x}^{
2}}}B_{{0}} \left( x \right) -f_{{B}}=0\cr
v{\frac {\rm d}{{\rm d}x}}C_{{0}} \left( x \right) &-  D  
{\frac {{\rm d}^{2}}{{\rm d}{x}^{2}}}C_{{0}} \left( x \right) -f_{{C}} = 0\cr
\alpha\,A_{{1}} \left( x \right) &-B_{{0}} \left( x \right) A_{{0}}
 \left( x \right) +v{\frac {\rm d}{{\rm d}x}}A_{{1}} \left( x \right) 
-  D  {\frac {{\rm d}^{2}}{{\rm d}{x}^{2}}}A_{{1}}
 \left( x \right)=0\cr 
\alpha\,B_{{1}} \left( x \right) &-B_{{0}} \left( x \right) A_{{0}}
 \left( x \right) +v{\frac {\rm d}{{\rm d}x}}B_{{1}} \left( x \right) 
-  D  {\frac {{\rm d}^{2}}{{\rm d}{x}^{2}}}B_{{1}}
 \left( x \right)=0\cr 
B_{{0}} \left( x \right) A_{{0}} \left( x \right) &+v{\frac {\rm d}{
{\rm d}x}}C_{{1}} \left( x \right) - D  {\frac {{\rm d}^
{2}}{{\rm d}{x}^{2}}}C_{{1}} \left( x \right) =0\cr 
\alpha\,A_{{2}} \left( x \right) &-B_{{0}} \left( x \right) A_{{1}}
 \left( x \right) -B_{{1}} \left( x \right) A_{{0}} \left( x \right) +
v{\frac {\rm d}{{\rm d}x}}A_{{2}} \left( x \right) -  D  
{\frac {{\rm d}^{2}}{{\rm d}{x}^{2}}}A_{{2}} \left( x \right)=0 \cr
\alpha\,B_{{2}} \left( x \right) &-B_{{1}} \left( x \right) A_{{0}}
 \left( x \right) -B_{{0}} \left( x \right) A_{{1}} \left( x \right) +
v{\frac {\rm d}{{\rm d}x}}B_{{2}} \left( x \right) -  D  
{\frac {{\rm d}^{2}}{{\rm d}{x}^{2}}}B_{{2}} \left( x \right) =0\cr
B_{{1}} \left( x \right) A_{{0}} \left( x \right) &+B_{{0}} \left( x
 \right) A_{{1}} \left( x \right) +v{\frac {\rm d}{{\rm d}x}}C_{{2}}
 \left( x \right) -  D  {\frac {{\rm d}^{2}}{{\rm d}{x}^{
2}}}C_{{2}} \left( x \right)=0 
\cr}$$
etc.  Each of these equations involves only one differentiated variable, occurring linearly and with constant coefficients, the other $A_i$, $B_i$ and $C_i$ occurring in the "inhomogeneous terms" being already known from the previous equations.  Thus they may be solved one by one.
A: There are a few things to notice about your system:


*

*By shifting $c_{A/B/C}$, i.e. introducing $\hat{c}_A = c_A + \delta_A$ and similarly for $c_{B,C}$, you can get rid of the constant source terms $f_{A/B/C}$.

*The first two equations are decoupled from the third one, so you can first focus on solving $c_{A}$ and $c_B$, and use that as input in the third equation.

*You can remove all the nonlinear terms by a change of variables. Consider $u = c_A - c_B$; when you substract the second model equation from the first, you see that the resulting equation for $c_A - c_B = u$ is linear. The same holds for $v = c_A + c_C$ and for $w = c_B + c_C$. Using these new variables, you will end up with three decoupled, linear equations for $u,v,w$. Moreover, the equations for $u,v,w$ are identical up to the source term $f$. Since there is a linear relation between $u,v,w$ (in this case, $u = v - w$), this procedure does not solve the system completely. However, it will get you started.

