% Lösning till uppgift ENM 7.13 -- Euler framåt för system
% Joar Bagge 2016-11-11

clear

h = 0.05; % steglängd
w0 = [1 0 0 1.4]'; % initialvärden (kolonnvektor)
t0 = 0; tend = 10; % tidsintervall

% Högerled till ekvationen w' = f(w)
% Obs: f är en funktion från R^4 till R^4.
f = @(w) [w(2);
          -w(1)/(w(1)^2+w(3)^2)^(3/2);
          w(4);
          -w(3)/(w(1)^2+w(3)^2)^(3/2)];

N = (tend-t0) / h; % antal steg
t = t0:h:tend; % vektor med alla tidpunkter
w = zeros(length(w0),length(t)); % lösningsmatris
w(:,1) = w0; % sätt initialvärdet!
% Obs: varje kolonn i matrisen w motsvarar lösningen vid en viss tidpunkt.

% Euler framåt
for n=1:N
    w(:,n+1) = w(:,n) + h*f(w(:,n));
end

% Plotta planetens bana, dvs (x(t), y(t))
plot(w(1,:), w(3,:))
hold on
plot(0,0,'*') % solens position
hold off
axis equal
