% ENM 7.15
% Denna kod skrevs efter övningen 2017-10-11 i SF1514/18/19

% 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 hjälp av inskjutningsmetoden. Detta betyder
% att vi i första steget vill bortse från villkoret vid x=1 och
% lösa diffekvationen som ett begynnelsevärdesproblem, med
% villkor vid x=0. Vi har bara ett villkor vid x=0 så vi måste
% hitta på ett annat.

clear

%%
% Första steget -- lös begynnelsevärdesproblem.
%
% Skriv om diffekvationen som ett system av första ordningen. Vi
% inför variablerna u1=y och u2=y' och får då
%
%   u1' = u2
%   u2' = g(u1) - (2/x)*u2    så länge x inte är noll.
%
% Om x=0 använder vi l'Hopitals regel för att få ekvationen
%
%   u2' = g(u1)/3.
%
% Detta är en ekvation på formen u' = f(x,u) där u=[u1;u2].
% Vi behöver begynnelsevillkor och vet att u2(0)=y'(0)=0, men
% saknar villkor för u1(0). Vi hittar på u1(0)=a. Alltså blir
% u(0) = [a; 0]. Vi provar olika värden på a och ser vad vi får
% för resultat:

clf
for a = [0, 0.5, 0.75, 1]

  u0 = [a, 0]; % begynnelsevillkor
  [X,U] = ode45(@f, [0, 1], u0); % lös u'=f(x,u) med ode45 på
                                 % intervallet 0 <= x <= 1,
                                 % med begynnelsevillkor u0.
  % Funktionen f som används i anropet till ode45 definieras
  % längre ner i den här filen.

  plot(X, U(:,1)), hold on % rita ut lösningen

end
disp('Titta på figuren som dök upp');
% Tydligen ger a=0.75 ett värde på y(1) som är ganska nära 1.
disp('Tryck på en tangent för att fortsätta');
pause;

%%
% Andra steget -- bestäm vilket a som behövs noggrannt
%
% Vi kan se processen att lösa begynnelsevärdesproblemet som en funktion F
% som tar a som indata och returnerar värdet på lösningen y i punkten x=1.
% Alltså F(a) = y(1). Om vi vill hitta det a som gör att y(1)=1 så ska vi
% lösa ekvationen F(a)=1. Det är samma sak som att hitta ett nollställe
% till funktionen F(a)-1, som vi kallar skjut(a) nedan: skjut(a)=F(a)-1.
%
% Vi kan använda valfri metod för att hitta nollstället till funktionen
% skjut, t.ex. sekantmetoden, intervallhalvering eller Matlabs inbyggda
% funktion fzero. (Newtons metod passar inte så bra eftersom vi inte kan
% derivera funktionen skjut för hand.) Vi använder fzero nedanför.

a0 = fzero(@skjut, 0.75) % hitta nollstället till skjut.
                         % startgissningen är 0.75 baserat på första steget.
% Funktionen skjut definieras längre ner i den här filen.

% Upprepa till slut lösningen med det rätta värdet på a så vi får hela
% lösningskurvan och även kan räkna ut eta.
u0 = [a0, 0];
[X,U] = ode45(@f, [0, 1], u0);
eta = 3*U(end,2) % 3*y'(1)


%%
% Funktion för första steget: u' = f(x,u)
function dudx = f(x, u)
  g = @(y) y.*exp(4*(1-y)./(1+0.2*(1-y)));
  if x == 0
    dudx = [u(2); g(u(1))/3];
  else
    dudx = [u(2); g(u(1)) - 2*u(2)/x];
  end
end

%%
% Funktion för andra steget. 'skjut' tar a som indata, löser
% begynnelsevärdesproblemet och returnerar skillnaden mellan lösningens
% värde i x=1 och det värde vi vill ha, nämligen 1.
function skillnad = skjut(a)
  u0 = [a, 0];
  [X,U] = ode45(@f, [0, 1], u0);
  plot(X, U(:,1)), hold on
  skillnad = U(end,1) - 1; % skillnad mellan y(1) och vad det borde vara
end
