% M-File grxxm044.m
% Euler-Verfahren / Heun-Verfahren
% y' = a y sin(at) , t > 0,   y(0) = y_0

clear;
figure('Position',[50 50 1000 800],'color',[1 1 0.77]);
%% 'Position' in Octave nicht definiert

% Differentialgleichung
f=inline('a*y*sin(a*t)','t','y','a');
% Loesung
l=inline('y0*exp(1-cos(a*t))','t','y0','a');
% Parameter
a=4.1;
y0=3.5;
T=2*pi;

% Schrittweite fuer Euler- und Heun-Verfahren
n=400;
h=T/n;
hh=num2str(h);

% Loesungsplot auf [0,2*pi]
subplot(3,1,1);
tt=0:T/1000:T;
yy=l(tt,y0,a);
plot(tt,yy,'b');
legend('Loesung');
aa=num2str(a);
y00=num2str(y0);
Ltext=['Loesung des AWP:    dy/dt = a y sin(at) ,   y_0 = ' y00 ',    (a = ' aa ')'];
title(Ltext,'FontSize',15);
xlabel('t - Achse');

% Euler-Verfahren 
subplot(3,1,2);
E(1,1:4)=[0 y0 y0 0];
for k=1:n
    E(k+1,1)=k*h;
    E(k+1,2)=E(k,2)+h*f(E(k,1),E(k,2),a);
    E(k+1,3)=l(E(k+1,1),y0,a);
    E(k+1,4)=E(k+1,3)-E(k+1,2);
end;
plot(E(:,1),E(:,2),'m',E(:,1),E(:,3),'b');
legend('Euler','Loesung');
Etext=['Euler-Verfahren:    h = ' hh ];
title(Etext,'FontSize',15);
xlabel('t - Achse');

% Heun-Verfahren
subplot(3,1,3);
H(1,1:4)=[0 y0 y0 0];
for k=1:n
    H(k+1,1)=k*h;
         k1=f(H(k,1),H(k,2),a);
         k2=f(H(k+1,1),H(k,2)+h*k1,a);
    H(k+1,2)=H(k,2)+h*(k1+k2)/2;
    H(k+1,3)=l(H(k+1,1),y0,a);
    H(k+1,4)=H(k+1,3)-H(k+1,2);
end;
plot(H(:,1),H(:,2),'r',H(:,1),H(:,3),'b');
legend('Heun','Loesung');
Htext=['Heun-Verfahren:    h = ' hh ];
title(Htext,'FontSize',15);
xlabel('t - Achse');

% Tabelle 
disp('---------------------------------------------------------------------');
disp('          Loesung        Euler-Verfahren          Heun-Verfahren ');
disp('---------------------------------------------------------------------');
disp('  t_n      y(t_n)        E_n      E_n-y(t_n)       H_n     H_n-y(t_n) ');
disp('---------------------------------------------------------------------');
% Ausgabe der ersten s Werte 
% Ausgabe aller Werte wenn 100 durch n ersetzt wird
for s=1:50
    fprintf('%2.5f %10.6f %12.6f %10.6f %13.6f %10.6f',...
        E(s,1),E(s,3),E(s,2),E(s,4),H(s,2),H(s,4));
    fprintf('\n');
end;

% Schliessen
%% unter Octave: nachfolgende Befehle kommentieren!
uicontrol(gcf,...
         'Style','Push',... 
         'Position',[440 10 120 40],...
         'String','Schliessen',...
         'Backgroundcolor',[0.9 0.93 0.31],...
         'callback','close');
