%% Valeur maximale de z et intervalle de valeurs de z
zmax= 10; zint= [0 zmax];

%% Resolution et traces pour divers parametres m de Falkner-Skan
figure; hold on; calcsol(0, 0.3, zint, 'k');
% calcsol(0.1, ..., zint, 'r');
% calcsol(0.2, ..., zint, 'm');

%% Fonction principale calcsol, alphaest est une estimation de alpha optimal
function calcsol(m, alphaest, zint,  color)
  disp(['m= ',num2str(m)]);

  % Calcul de alpha avec l'objectif que tir = 0
  % Le 1er argument est la fonction a annuler,
  % le 2d donne l'intervalle dans lequel chercher alpha
  alpha= fzero(@(x) tir(x , zint, m) , [0.8*alphaest 1.2*alphaest])

  % Resolution de l'equation differentielle une fois alpha connu
  [z,F]= ode45(@(z,F) G(z,F,m) , zint , [ 0 ; 0 ; alpha ]);

  % Recuperation des valeurs tabulees de la fonction et de sa derivee
  f= F(:,1); fprime= F(:,2);

  % Graphe de f'(z)
  plot(z, fprime, color, 'LineWidth',2);
  xlabel("z"); ylabel("fprime"); title(['m = ',num2str(m)]);
end

%% Fonction tir qui donne la valeur de f'-1 en zmax en fonction de alpha
function res= tir(alpha , zint, m)
  disp(['alpha= ',num2str(alpha)]);
  [z,F]= ode45(@(z,F) G(z,F,m) , zint , [ 0 ; 0 ; alpha ]);
  res= F(end,2) - 1;
  disp(['fprime(zmax) - 1 = ',num2str(res)]);
end

%% Fonction G qui definit l'equation differentielle
function res= G(z, F, m)
  res= [F(2) ; F(3) ; ... ];
end
