% ENM 7.11
% Denna kod skrevs till övningen 2017-11-01 i SF1545

clear

% Vi ska här lösa en till andra ordningens diffekvation, nu med
% ode45. Vi måste skriva om den som ett system av första
% ordningens diffekvationer först. De två variablerna blir höjden
% h och hastigheten dh/dt.

H = 2000; % starthöjd

% Tidintervall (start och slut)
t_start = 0;
t_end = 60*3;

y0 = [H; 0]; % initialvillkor

[t,y] = ode45(@f, [t_start t_end], y0);
% ode45 returnerar en vektor med alla tidpunkter och en matris
% med lösningen, där varje kolumn motsvarar en variabel och varje
% rad en tidpunkt (tvärtom som vi gjorde i K6.5).

% Vi måste skriva @f för att Matlab ska förstå att vi skickar in
% själva funktionen f till ode45, utan snabel-a:et skulle Matlab
% försöka evaluera funktionen redan i anropet (det blir ett
% felmeddelande).

plot(t, y(:,1)) % rita första variabeln, alltså höjden
grid on


% Funktion för högerledet i diffekvationen y'=f(t,y)
function dydt = f(t, y)
  % Konstanter givna i uppgiften
  m = 75;
  g = 9.81;
  vf = 50;
  vs = 5;
  % Från sluthastigheterna vf och vs kan man räkna ut
  % luftmotståndskoefficienten k genom att sätta in i den
  % ursprungliga diffekvationen (se uppgiften) att
  % andraderivatan d2h/dt2 är noll och hastigheten dh/dt är
  % sluthastigheten. Lös sedan ut k.
  if t < 30
    k = m*g/vf^2;
  else
    k = m*g/vs^2;
  end
  dydt = [y(2); (k/m)*y(2)^2 - g];
end
