% ENM 7.15
% Denna kod skrevs på övningen 2017-11-09 i SF1545

% Problemet vi vill lösa är att bestämma y(x) så att
%
%   y'' + (2/x)*y' = g(y)     för 0 <= x <= 1
%   där g(y) = y*exp(4*(1-y)/(1+0.2*(1-y)))
%
%   y'(0) = 0  och  y(1) = 1.
%
% Vi gör detta med inskjutningsmetoden. Se slides för en
% sammanfattning av hur den funkar.

clear

% Första steget är att skriva om problemet som ett
% begynnelsevärdesproblem. Vi saknar ett begynnelsevärde och
% måste därför gissa. I detta fall gissar vi att
% y(0) = a = 0.5.

a = 0.5; % gissning
u0 = [a, 0]; % startvärde
[t,u] = ode45(@fun, [0 1], u0) % lös med ODE45
figure(1)
plot(t,u)

% Vi kollar hur väl uppfyllt y(1) = 1 är. Inte så bra i detta
% fall. Vi kan testa att ändra a manuellt, men mer strukturerat
% är att se det som ett ekvationslösningsproblem. Hitta på
% funktionen F(a) så att F(a)=0 om slutvärdet är uppfyllt.

a0 = fzero(@F, 0.7) % hitta nollställe till F med fzero (startgissning 0.7)
F(a0) % är funktionen noll?
u0 = [a0, 0];
[t,u] = ode45(@fun, [0 1], u0);
figure(2)
plot(t,u) % kolla att villkoret är uppfyllt nu

% BONUS: Hur kan man hitta en startgissning till fzero förutom
% att gissa? Jo, man kan plotta funktionen F(a).
avec = linspace(0, 2, 100);
% Det går inte att skriva Fvec = F(avec) eftersom funktionen F
% förväntar sig att indata är en skalär. Man kan istället skriva
Fvec = arrayfun(@F, avec);
% vilket är samma sak som att skriva en for-loop och applicera F
% på ett a-värde i taget (funktionen arrayfun gör alltså precis
% det).
figure(3)
plot(avec, Fvec), xlabel('a'), ylabel('F(a)')
grid on
% I figuren ser vi att a=0.7 är ett bra värde.


% Funktion för högerledet till begynnelsevärdesproblemet.
function dudx = fun(x,u)
    g = @(y) y*exp(4*(1-y)/(1+0.2*(1-y)));
    if x > 0
        dudx = [u(2);
                g(u(1)) - (2/x)*u(2)];
    elseif x == 0
        dudx = [u(2);
                g(u(1))/3];
    end
end

% Funktion som visar hur väl slutvärdet är uppfyllt.
function slutfel = F(a)
    u0 = [a, 0];
    [t,u] = ode45(@fun, [0 1], u0);
    y1 = u(end,1); % motsvarar y(1)
    slutfel = y1 - 1;
end
