% freie gedämpfte/ungedämpfte Schwingung

clear all
hold off

% Richtungsfeld:
% Punkte, wo ein Richtungsvektor angebracht werden
% soll:

  xcord = -1:.1:1;
  ycord = -1:.1:1;

  [xx,yy] = meshgrid(xcord,ycord);

% Bestimmung der Richtungsvektoren in den
% entsprechenden Punkten: (geg. Dgl.)

  delta = 0.0;
  xxp = yy;
  yyp = - delta.*yy-xx;

  quiver(xx,yy,xxp,yyp)
  grid
  pause

  hold on

% Lösung der Dgl.
% Anfangszeit: 
  t0= 0;
% Endzeit: 
  tfinal = 7;
  
  tspan = [t0 tfinal];
% Definition der AB:

  x0 = [0 1]';

% eine Dgl 2-ter Ordnung muss auf ein System von 
% Dgl 1-ter Ordnung umgeschrieben werden


[t, x] = ode23('dgl3', tspan, x0);

% dabei ist folgende Interpretation gültig:
% x1=x;
% x2=xp;


plot(x(:,1),x(:,2), 'r',  'LineWidth', 2)

pause
hold on

axis equal
grid
title(['Phasenporträt einer freien gedämpften Schwingung', ' delta = ', num2str(delta)])
xlabel('x (=Ort)')
ylabel('xp (=Geschwindigkeit)')
grid

ta = 0;
te = tfinal;
h = 0.1;  % h = 0.5;
ya = [0 1]';
[x, y] = euler_kepler('dgl3',ta,te,h,ya);
plot(y(1,:), y(2,:), 'k.-', 'LineWidth', 2)

hold on
pause

ta = 0;
te = tfinal;
h = 0.1;  % h = 0.5;
ya = [0 1]';
[x, y] = euler_implizit('dgl3i', ta, te, h, ya);
plot(y(1,:), y(2,:), 'm.-', 'LineWidth', 2)

hold on

pause

ta = 0;
te = tfinal;
h = 0.1;  % h = 0.5;
ya = [0 1]';
[x, y] = euler_implizit('dgl3imMP', ta, te, h, ya);
plot(y(1,:), y(2,:), 'g.-', 'LineWidth', 2)

legend('VF', 'exakt', 'Euler explizit', 'Euler implizit', 'imMP')

hold off