% gestörter Harmonischer Oszillator

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

hold off
clear all, clc

  xcord = -3:.5:3;
  ycord = -3:.5:3;
  
%   xcord = -10.5:1:10.5;
%   ycord = -5.5:1:5.5;

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

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

  eps = .1;
  xxp = yy;
  yyp = + eps.*xx.^3-xx;

  quiver(xx,yy,xxp,yyp)
  grid

  hold on

% Lösung der Dgl.
% Anfangszeit: 
  t0= 0;
% Endzeit: 
  tfinal=20;

% Definition der AB:

  x0 = [0 1]';

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

tspan = [t0, tfinal];

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

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


plot(x(:,1),x(:,2), 'Color', [0 0.4*norm(x0) .1]);
title(['gestörter harmonischer Oszillator: eps = ', num2str(eps)])
xlabel('x (=Ort)')
ylabel('xp (=Geschwindigkeit)')
grid
hold on

x0 = [0 2]';
[t, x] = ode23('dgl3e', tspan, x0);
plot(x(:,1),x(:,2), 'Color', [0 0.4*norm(x0) .1]);
hold on

x0 = [0 0.5]';
[t, x] = ode23('dgl3e', tspan, x0);
plot(x(:,1),x(:,2), 'Color', [0 0.4*norm(x0) .1]);
hold on

x0 = [0 1.5]';
[t, x] = ode23('dgl3e', tspan, x0);
plot(x(:,1),x(:,2), 'Color', [0 0.4*norm(x0) .1]);
grid on


