Symmetric coupled second order differential equations

Question

I have two coupled second order differential equations that seem very similar to this related problem (A heat transfer textbook, https://ahtt.mit.edu/),

$$ (b-(a+x))y_1''-y_1'+c(y_1-y_2) = 0\\ (b+(a+x))y_2''+y_2'+c(y_2-y_1) = 0, $$ that I need to solve on the domain $-a\leq x\leq a$. $y_1$ and $y_2$ are dimensionless steady state temperature profiles that are only a function of $x$. I think that the profiles $y_1(x) = -y_2(-x)$ because of the symmetry in the problem if the boundaries are $y_1(-a)=-y_2(a)=T_\infty$ and $y_1'(-a)=y_2'(a)$. This would lead me to believe that I only need to solve one of the two equations. Unfortunately, quoting from the textbook, "this equation is difficult to solve because it has a variable coefficient" and I've just made it a bit more complex. Could you tell me

• if it's possible to solve this system of equations analytically? I have tried solving by replacing $y_2 \rightarrow -y_1(-x)$ and guessing some combinations of exponential functions for $y_1$ such that I could get rid of the $c$-term. The Bessel function in the textbook suggests it can't be that simple. How would I start solving the problem? Or;
• if it's possible to solve the system numerically? I am familiar with time propagating ODE solvers in Matlab, but here I have a steady state problem with boundary conditions rather than initial conditions. Any software you would recommend?

Answer 1

If the boundary conditions are $y_1(-a)=-y_2(a)=T_{\infty}$ and $y_1'(-a)=y_2'(a)$, one can solve the system of ODEs with the ansatz $y_1(x)=A_1x+B_1, y_2(x)=A_2x+B_2$. Indeed,

1. The boundary condition $y_1'(-a)=y_2'(a)$ is satisfied if $$ A_1=A_2=A. \tag{1} $$
2. Plugging the ansatz into the ODEs we obtain \begin{align} -A+c(B_1-B_2)&=0, \tag{2a} \\ A+c(B_2-B_1)&=0. \tag{2b} \end{align} Both equations are satisfied if $$ B_1-B_2=\frac{A}{c}\quad(c\neq 0). \tag{3} $$
3. The boundary conditions $y_1(-a)=-y_2(a)=T_{\infty}$ imply \begin{align} -Aa+B_1&=T_{\infty}, \tag{4a} \\ -Aa-B_2&=T_{\infty}. \tag{4b} \end{align} Subtracting $(4b)$ from $(4a)$ we obtain $$ B_1+B_2=0. \tag{5} $$
4. It follows from $(3)$ and $(5)$ that $$ B_1=-B_2=\frac{A}{2c}. \tag{6} $$
5. Finally, plugging $(6)$ into $(4a)$ we obtain $$ A=\frac{T_{\infty}}{-a+\frac{1}{2c}}\quad(2ac\neq 1). \tag{7} $$
6. In conclusion, we've shown that the solution to the given system of differential equations satisfying the boundary conditions $y_1(-a)=-y_2(a)=T_{\infty}$ and $y_1'(-a)=y_2'(a)$ is given by \begin{align} y_1(x)&=\frac{T_{\infty}}{-a+\frac{1}{2c}}\left(x+\frac{1}{2c}\right)= T_{\infty}\left(\frac{2cx+1}{1-2ca}\right), \tag{8a} \\ y_2(x)&=T_{\infty}\left(\frac{2cx-1}{1-2ca} \tag{8b}\right), \end{align} provided $2ca\neq 1$ (the condition $c\neq 0$ in $(3)$ can be relaxed in the final solution).