% M-File grxxm054.m
% Euler-Runge-Kutta-Verfahren
% Feder-Masse-system

figure('Position',[50 50 1000 800],'color',[1 1 0.77]);
% 'Position' unter Octave nicht definiert (Warnung!).
clf;
clear;

% Testbeispiel
z1=inline('sqrt(3)/6*sin(t*sqrt(3)/sqrt(m))+sin(t/sqrt(m))/2','t','m');
z2=inline('sin(t/sqrt(m))/2-sqrt(3)/6*sin(t*sqrt(3)/sqrt(m))','t','m');

m=2.5;
T=6*pi;

% ODE Systemmatrix
A=[0 0 1 0;0 0 0 1;-2/m 1/m 0 0;1/m -2/m 0 0];

% Plot Testloesung
tt=0:T/1000:T;
z1p=z1(tt,m);
z2p=z2(tt,m);

subplot(3,1,1);
plot(tt,z1p,'b',tt,z2p,'g',[0,T],[0,0],'k');
legend('z_1(t)','z_2(t)');
xlabel('t - Achse');
title('Loesung des Testbeispiels','FontSize',15);

% Runge-Kutta-Verfahren
n=40;   %  Schrittweitenwahl
h=T/n;
hh=num2str(h);
R(1,1:5)=[0 0 0 1/sqrt(m) 0];  % Anfangsbedingungen
for k=1:n
    R(k+1,1)=k*h;
    w=R(k,2:5)';
    k1=A*w;
    k2=A*(w+h/2*k1);
    k3=A*(w+h/2*k2);
    k4=A*(w+h*k3);
    u=w+(k1+2*k2+2*k3+k4)*h/6;
    R(k+1,2:5)=u';
end;

% Plot der R-K-Naeherungen
subplot(3,1,2);
plot(R(:,1),R(:,2),'r',tt,z1p,'b',[0,T],[0,0],'k');
legend('z_1 - R-K-Verfahren','z_1(t)');
xlabel('t - Achse');
z1Text=['Runge-Kutta-Verfahren:   z_1(t)  mit   h = ' hh];
title(z1Text,'FontSize',15);

subplot(3,1,3);
plot(R(:,1),R(:,3),'m',tt,z2p,'g',[0,T],[0,0],'k');
legend('z_2 - R-K-Verfahren','z_2(t)');
xlabel('t - Achse');
z1Text=['Runge-Kutta-Verfahren:   z_2(t)  mit   h = ' hh];
title(z1Text,'FontSize',15);


% Schliessen
% (unter OCTAVE arbeitet der folgende Befehl des GUI nicht)
% Octave: Man ersetze die naechsten 5 Zeilen durch Kommentare
uicontrol(gcf,...
         'Style','Push',... 
         'Position',[440 10 120 40],...
         'String','Schliessen',...
         'Backgroundcolor',[0.9 0.93 0.31],...
         'callback','close');
