How can I solve this 2nd-order, non-linear, stiff ODE numerically using Matlab?

Question

I need to find a numerical solution for $T(x)$ to

\begin{equation} aT'' - bT^4 + c = 0, \ \ \ T(0) = T_0, \ \ \ T'(L) = 0 \end{equation} where $a$, $b$, $c$, $T_0$ and $L$ are positive constants. This is a heat conduction and emission problem with boundary conditions, and FYI we know that $T'(0)<0$.

The problem is stiff, and Mathematica is failing me.

What Matlab solver do I use? How do I find an equivalent system of first-order equations?

Answer 1

An example:

b_over_a=1;
c_over_a=1;
T0=2;
L=1;
odefun=@(x,y) [y(2);b_over_a*y(1).^4-c_over_a];
bcfun=@(ya,yb) [ya(1)-T0;yb(2)];
solinit=bvpinit(linspace(0,L,1e2),[T0;0]);
sol=bvp4c(odefun,bcfun,solinit);
plot(sol.x,sol.y(1,:))

This plots the solution to your problem with $b/a=1,c/a=1,T_0=2,L=1$. It uses a constant initial guess i.e. $T(x)=2,T'(x)=0$, and uses an initial mesh with 100 evenly spaced points. It gives Matlab the option to extend this mesh, but in this particular case it decides not to do so.

Your problem might be stiffer than this one depending on the exact parameter values you are using, but the particular parameters that I chose appear to not be stiff at all.

For more details see http://www.mathworks.com/help/matlab/ref/bvp4c.html

Answer 2

If this helps,

$$aT'T'' +(c - bT^4)T' = 0$$

then

$$\frac a2T'^2+(cT-\frac b5T^5)=C.$$