% Der harmonische Oszillator - symplektisches Eulerverfahren

T  = 100;         % Stop-Zeitpunkt
dt = 0.1;         % Zeitschritt
N  = floor(T/dt); % Dann brauchen wir N Schritte

x = zeros(2,N);   % Fuer bessere Effizienz preallokiere Loesungsvektor
x(:,1) = [0; 1];  % Anfangsdaten

M1 = [1 0; dt 1];
M2 = [1 dt; 0 1];
M  = inv(M1) * M2;

for t=1:N-1
x(:,t+1) = M * x(:,t);
end

plot(x(1,:),x(2,:)), grid, replot 
