function tpbvp

global c h Tinf

% Constants
k = 42;
b = 0.5;
d = 0.2;
Tinf = 500;
l = 2;
epsilon = 0.1;
sigma = 5.7e-8;
h = 0.5;
Tr = 675;

% Computed constants
A = b*d;
P = 2*(b+d);
c = (sigma*epsilon*P)/(k*A);

% Boundary conditions
T0 = 1000;
T4 = 350;

% Part a)
% Finite difference approach with nonlinear equations and rootfinding

% Initial guess
T = zeros(3,1);
Tinit = zeros(3,1);
Tinit(1) = Tr;
Tinit(2) = Tr;
Tinit(3) = Tr;

% Nonlinear root finding
Tfd_nlin = fsolve(@finitediff,Tinit,optimset('fsolve'))

% T =
%   739.9453
%   592.5973
%   474.1394

% Part b)
% Finite difference approach with linearized system of equations

diag = 2+4*c*Tr^3*h^2;
rhs = c*(3*Tr^4+Tinf^4)*h^2;

Amat = [diag -1 0; -1 diag -1; 0 -1 diag];
bvec = [T0+rhs; rhs; T4+rhs];

Tfd_lin = Amat\bvec

% Part c)
% Shooting method with rootfinding

% Initial guess
Sinit = -800;

% Nonlinear root finding
S = fzero(@shoot,Sinit)
yinit = [S; T0];
[x,y] = numinteg(yinit);
Tsh_nlin = y(2:4,2)

%------------------------------------------------%

function zero = finitediff(T)

T1 = T(1);
T2 = T(2);
T3 = T(3);

zero(1) = -2*T1-0.0475e-8*T1^4+T2+1029.6875;
zero(2) = T1-2*T2-0.0475e-8*T2^4+T3+29.6875;
zero(3) = T2-2*T3-0.0475e-8*T3^4+379.6875;

%------------------------------------------------%

function zero = shoot(Sinit)

Tinit = 1000;
yinit = [Sinit; Tinit];

[x,y] = numinteg(yinit);

zero = y(5,2)-350;

%------------------------------------------------%

function [x,y] = numinteg(yinit)

xspan = linspace(0, 2, 5);
[x,y] = ode45(@eqn,xspan,yinit);

%------------------------------------------------%

function dydx = eqn(x,y)

global c h Tinf

S = y(1);
T = y(2);

dSdx = c*(T^4-Tinf^4);
dTdx = S;

dydx = [dSdx; dTdx];