function [x,u,q] = varheatSHOOT()
% Set up and solve, by nonlinear shooting, a "serious" linear
% two-point BVP for equilibrium temperature distribution in a rod.
% The terms correspond to variable conductivity, constant
% chemical-reaction-created heat, and variable externally-introduced heat.
% Uses bisection to solve target boundary equation u(3) = 0.
% Example:
% >> [x,u,q] = varheatSHOOT();
L = 3;
k = @(x) 0.5 * atan((x-1.0) * 20.0) + 1.0;
s = @(x) exp(-(x-2.0).^2);
r0 = 0.5;
% ODE Y' = G(x,Y) is described by this right-hand side
G = @(x,Y) [- Y(2) / k(x); % Y(1) = u
r0 * Y(1) + s(x)]; % Y(2) = q
% bracket unknown u(0)
a = -10.0; % produces u(3) which is too high
b = 0.0; % ... u(3) which is too low
% unsophisticated bisection method:
N = 100;
for n = 1:N
fprintf('.')
c = (a+b)/2;
% evaluate F(c) = (estimate of u(3) based on u(0)=c)
[xout,Y] = ode45(G,[0.0 3.0],[c; 0.0]);
F = Y(end,1);
if abs(F) < 1e-12
fprintf('\nbisection done; F(A)=%2e using A=%5f\n', F, c)
break % we are done
elseif F >= 0.0
a = c;
else
b = c;
end
end
if n >= N, fprintf(' bisection: max iterations reached\n'), end
% redo to get final version on a grid for plot
x = 0:0.05:3.0;
[xout,Y] = ode45(G,x,[c; 0.0]);
u = Y(:,1)';
q = Y(:,2)';
figure(2), plot(x,k(x),'r',x,s(x),'b',xout,u,'g*',xout,q,'k')
grid on, xlabel x, legend('k(x)','s(x)','u(x)','q(x)')
title(sprintf('result of SHOOTING: u(0) = %.6f',u(1)))