% Der harmonische Oszillator - implizites 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

A = [0 1; -1 0];
M = eye(2) - dt * A;
IM = inv(M);

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

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