function ShootingMethod
% Simple Matlab program to demonstrated the shooting method
% with the heat transfer example problem used in lecture.
% The differential equation for this problem is
% d^2T/dx^2 + h(Ta - T) = 0
% with x = 0 to 10, Ta = 20, and boundary conditions
% T(0) = 40 and T(10) = 200.

% First convert the second order ODE into two first order ODEs
% dT/dx = z (first derivative of original variable = new variable)
% dz/dx = d^2T/dx^2 = - h(Ta - T)
% The state vector y is then
% y = [T z]
% and the state derivative vector is
% dy/dx = [dT/dx dz/dx]

% For the shooting method, guess dT/dx(0), see if it gives the
% correct T(10) value, and then iterate until converged

% Approach 1: Iterate manually 5 times
T0 = 40;
T10 = 200;
npts = 101;
xspan = linspace(0,10,101);
    
for i = 1:5
    z0 = input('\nEnter a guess for dT/dx(0): ');
    y0 = [T0 z0]';
    [x,y] = ode45(@HeatTransfer,xspan,y0);

    T = y(:,1);
    fprintf('The final value of T is %5.2f\n', T(npts))
end

fprintf('The final guess for dT/dx(0) was %5.2f\n', z0)
plot(x,T)
xlabel('x')
ylabel('T')
pause

% Approach 2: Perform two iterations and use a linear curve fit
% First enter a guess for dT/dx(0) that will produce a value of T(10)
% that is too low (use dT/dx(0) = 10)
zlow = input('\nEnter a guess for dT/dx(0) such that T(10) is too low: ');
y0 = [T0 zlow]';
[x,y] = ode45(@HeatTransfer,xspan,y0);

Tlow = y(npts,1);
fprintf('The corresponding value of T(10) is %5.2f\n', Tlow)

% Next enter a guess for dT/dx(0) that will produce a value of T(10)
% that is too high (use dT/dx(0) = 20)
zhigh = input('\nEnter a guess for dT/dx(0) such that T(10) is too high: ');
y0 = [T0 zhigh]';
[x,y] = ode45(@HeatTransfer,xspan,y0);

Thigh = y(npts,1);
fprintf('The corresponding value of T(10) is %5.2f\n', Thigh)

% Now use a linear curve fit to find the value of dT/dx(0) that will make
% T(10) = 200. A linear curve fit works here since the ODE is linear.
slope = (zhigh-zlow)/(Thigh-Tlow);
z0 = zlow+slope*(T10-Tlow);
y0 = [T0 z0]';
[x,y] = ode45(@HeatTransfer,xspan,y0);
T = y(:,1);
fprintf('The final value of T is %5.2f\n', T(npts))
fprintf('The final value of dT/dx(0) is %5.2f\n', z0)

function dy = HeatTransfer(x,y)

% Define problem constants
h = 0.01;
Ta = 20;

% Unpack state vector
T = y(1,1);
z = y(2,1);

% Calculate state derivative vectory
dTdx = z;
dzdx = -h*(Ta-T);

% Save state derivative vector
dy(1,1) = dTdx;
dy(2,1) = dzdx;