function Hausaufgabe04
  h = [2, 1, 0.1, 1e-2, 1e-3]; % Schrittweiten
  x_exakt = 0:0.001:40;        % exakte Lösung erstellen
  y_exakt = exakt(x_exakt);
  
  fprintf('klassisches RKV\n');
  for i=1:5
      yrkv = klassRKV(0,40,h(i),[1;0],@(x,y) rhs(x,y)); % RK-Verfahren aufrufen
      
      %yeul = expleuler(0,40,h(i),[1;0],@(x,y) rhs(x,y)); % expliziter Euler zum Vergleich
     
      x = 0:h(i):40;
      %xeul = 0:h(i):40;
     
      figure(1);
      plot(x,yrkv(1,:),x_exakt,y_exakt);
      %figure(2);
      %plot(xeul,yeul(1,:),x_exakt,y_exakt);
      
      fprintf('Maximaler Fehler: %e\n', max(abs(yrkv(1,:)-exakt(x))));
      %fprintf('Maximaler Fehler euler: %e\n', max(abs(yeul(1,:)-exakt(xeul))));
      pause;
  end
end

function f = rhs(x,y) % rechte Seite nach Umformung in System erster Ordnung
  f = zeros(length(y),1);
  f(1) = y(2);
  f(2) = -0.12*y(2) - 2*y(1);
end

function f = exakt(x) % exakte Lösung
  f = exp(-0.06*x).*cos(sqrt(1.9964)*x)+.06/sqrt(1.9964)*exp(-0.06*x).*sin(sqrt(1.9964)*x);
end

function y = klassRKV(a,b,h,y0,f)
  x = a:h:b;
  n = length(x);
  m = length(y0); % Dimension des Systems bestimmen
  y = zeros(m,n); % Speicher bereitstellen
  y(:,1) = y0;    % Anfangswert setzen
  
  for i=1:n-1 % für jeden Zeitschritt: Zwischenwerte k_i berechnen, daraus dann neue Näherung y_i+1
      k1 = f(x(i),y(:,i));
      k2 = f(x(i)+1/2*h,y(:,i)+h*1/2*k1);
      k3 = f(x(i)+1/2*h,y(:,i)+h*1/2*k2);
      k4 = f(x(i)+h,y(:,i)+h*k3);
      y(:,i+1) = y(:,i) + h*(1/6*k1+1/3*k2+1/3*k3+1/6*k4);
  end
end

function y = expleuler(a,b,h,y0,f)
  x = a:h:b;
  n = length(x);
  m = length(y0);
  y = zeros(m,n);
  y(:,1) = y0;
  
  for i=1:n-1
      y(:,i+1) = y(:,i) + h*f(x(i),y(:,i));
  end
end
