%% uppgift 8

close all

g0 = 40;
gL = 60;
L = 1;
k = 2;

s = linspace(-10,110); % Guesses for unknown u' that will result in correct upper boundary value u(b,s)=gL. 
xspan = [0 1]; % Problem domain, NOT boundary conditions.

q = zeros(size(s));
for k = 1:length(s)
    [x, u] = ode45(@ode_system, xspan, [g0; s(k)]);
    q(k) = u(end,1) - gL; % Compare obtained boundary value with actual one.
end
plot(s,q)
hold on
title('Error in upper bound (q(s)) for various guesses for initial conditions (s) ')
xlabel('s')
ylabel('q(s)')
grid on

q_root = fzero(@(s)q_heat_transfer(s), 0) % Find initial value s that results in boundary satisfied.
plot(q_root, 0, 'ro')

figure 
[x u] = ode45(@ode_system, xspan, [g0; q_root]);
plot(x, u(:,1))
title('Heat transfer')
xlabel('x'); ylabel('u(x)')
hold on
%plot(x, u(:,2))
grid on
axis([0 1 40 65])

%

function ret = f_source(x)
    ret = 200*exp(-(x-1/2).^2);
end

function q = q_heat_transfer(s)
    [x u] = ode45(@ode_system, [0 1], [40; s]);
    q = u(end,1) - 60;
end

function uprim = ode_system(x,U)
    uprim = [U(2); 200*exp(-(x-1/2).^2)./(-2)];
end