2nd Order Non-Linear ODE in MATLAB Issues

[TurboRegal]

Hey everyone,
Having some trouble here using the solver we were supplied and modifying it to fit our problem...

I have a wire with a current flowing through it. I'm trying to find the temperature distribution wrt. position in the wire.

BV's are:
T(x=L/2) = 300K
dT/dx (x=0) = 0
(Apparently this is suppose to result in a symmetrical distribution around 0 from -L/2 -> L/2)The DE is as follows:
http://img165.imageshack.us/img165/3875/2ndorderodetm7.jpg

Which I've rearranged to:
d²y₂/dx² = phi1 * y₁ + phi2 * (y₁)⁴ - constants

Where y1 = T, y2 = dT/dx, and phi1, phi2, constants are just the collected constant terms...

Here is the code I'm using, in it c->T, r->x, sorry the terms were different and I haven't had time to change them yet...

%ThieleBVP.m - Example 4
clear
clc
global phi
global phi2
global constant

%Maximum number of secant iterations
maxsec = 20;

%Absolute error tolerance
tol = 1e-15;

%Arrays to store the guessed boundary condition, c(0), and the residual, R
c0 = zeros(maxsec, 1);
R = zeros(maxsec, 1);

%Set the phi is for T, phi2 is for T^4, constant is the rest of the
%collected terms
%Setting Constants...
kcon = 72
h = 2000
eps = 0.1
sig = 5.67 * 10^-8
Icur = 2
rowe = 32*10^-8
Tinf = 300
D = 80*10^-6
L = 4*10^-3

%Calculating Stuff...
RHS = -(Icur^2*rowe)/(kcon*(pi/4*D^2)^2)
LHS1 = -4*h/(kcon*D)
LHS2 = -4*eps*sig/(kcon*D)

%Finding the final constants...
phi = -LHS1
phi2 = -LHS2
constant = RHS + LHS1*Tinf + LHS2*Tinf

%Target boundary condition to shoot for, c(1)
c1_true = 300;

%First guess of the unknown boundary condition, c0 = c(0)
c0(1) = 280;

%Shooting loop
%rspan = [0 1]; % integration interval
rspan = [-L/2 L/2];

% first shot
i = 1;
yini = [c0(i) 0]; % initial conditions
[r,y] = ode45('ODEs',rspan,yini);
R(i) = y(length(r),1) - c1_true;

%second shot
i = 2;
c0(i) = 1.05*c0(i-1);
yini = [c0(i) 0]; % initial conditions
[r,y] = ode45('ODEs',rspan,yini);
R(i) = y(length(r),1) - c1_true;

while abs(R(i))>tol
i = i+1;
c0(i) = c0(i-1)-R(i-1)*(c0(i-1)-c0(i-2))/(R(i-1)-R(i-2));
yini = [c0(i) 0]; % initial conditions
[r,y] = ode45('ODEs',rspan,yini);
R(i) = y(length(r),1) - c1_true;
if i == maxsec
display('Maximum no. of iterations reached!');
break;
end
end

plot(r,y(:,1),'ro');

function dydr = ODEs(r,y)
global phi
global phi2
global constant
dydr = zeros(2,1);

dydr(1) = y(2);
dydr(2) = phi*y(1)+phi2*(y(1))^4-constant;

This just results in MATLAB spitting out jibberish.

At the moment, I need to figure out the following:
1) Does ode45 work for nonlinear ODE's?
2) How do I set the dT/dx (x=0) = 0 boundary value? I know this isn't working because if I set the phi2 power above to ^2 instead of ^4, I get:
http://img379.imageshack.us/img379/6339/80553833xt7.jpg
3) Why isn't MATLAB solving this properly?
4) Is the "solver" part of the 1st code even going to work for this type of problem. It was used at first for a 2nd order linear ODE, I'm assuming this is what is causing jibberish?

Thanks for the help...

 

[TurboRegal]

Fixing a few small errors in my stupid code...

%Finding the final constants...
phi = -LHS1
phi2 = -LHS2
constant = RHS + LHS1*Tinf + LHS2*Tinf^4

I think that's all I changed...

Nets me:
http://img292.imageshack.us/img292/453/18650296io7.jpg
http://g.imageshack.us/g.php?h=292&i=18650296io7.jpg


And gives me the error:
Warning: Failure at t=-1.512856e-003. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.469447e-018) at time t.
> In ode45 at 355
In ThieleBVP at 70

When I try just going from 0->L/2 and then making is symmetric myself it nets me the following:

This nets me (NO ERRORS):
http://img390.imageshack.us/img390/2945/99880853fh6.jpg
http://g.imageshack.us/g.php?h=390&i=99880853fh6.jpg

But it should be going up in the middle no?

 

[Jitse Niesen]

1) Does ode45 work for nonlinear ODE's?
2) How do I set the dT/dx (x=0) = 0 boundary value? I know this isn't working because if I set the phi2 power above to ^2 instead of ^4, I get:
http://img379.imageshack.us/img379/6339/80553833xt7.jpg

1. Yes

2. If you set TSPAN (the second argument of ode45) to [-L/2 L/2] then Matlab applies the initial conditions at -L/2 and integrates from x = -L/2 to x = L/2. Since you want to apply the initial conditions at x = 0, you have to set TSPAN to [0 L/2] (as you do below). Then you can compute the other half of the solution by symmetry, or by doing another integration but with TSPAN set to [0 -L/2].

When I try just going from 0->L/2 and then making is symmetric myself it nets me the following:

This nets me (NO ERRORS):
http://img390.imageshack.us/img390/2945/99880853fh6.jpg
http://g.imageshack.us/g.php?h=390&i=99880853fh6.jpg

But it should be going up in the middle no?

What do you mean with it should be going up in the middle? The solution looks fine to me, except that you get a temperature of -600K at x=0 which is physically impossible. Perhaps make sure that your constants are correct?