% Num. Math. - FS 2013 - Serie 4 Aufgabe 3 - Heun method

function heun
clear all
close all
%Anfangsbedingung
y0=[0.5;0;0;sqrt(3)];
%Endzeitpunkt
tf=500;
%Vektor mit versch. Schrittweiten fuer das Verfahren
%von Heun
h_vec=[2/10; 1/10; 1/100];
%Iteration ueber die Schrittweiten
for i=1:length(h_vec)
 h=h_vec(i); %Schrittweite
 N=tf/h;     %Anzahl Schritte
 y=zeros(length(y0),N+1); %Initializierung der Loesungs'matrix'
 y(:,1)=y0; %Initialisierung der Anf.bedingung
 %Iteration in der Zeit
 for k=1:N
  y(:,k+1)=...
 end
 %Plotten der ersten beiden Komponenten
 subplot(ceil(length(h_vec)/2)+1,2,i)
 plot(y(1,:),y(2,:))
 axis equal;
 axis([-3 1.5 -1.1 1.5])
 title(['Heun mit h=',num2str(h)]) 
end
 %Vergleich mit ode23
 [t,y]=... %!!!!!! Ergaenze hier den MATLAB ODE Solver !!!!!!
 subplot(ceil(length(h_vec)/2)+1,2,length(h_vec)+1)
 plot(y(:,1),y(:,2))
 axis equal;
 axis([-3 1.5 -1.1 1.5])
 title(['ode23']) 
 %Plotten der Schrittweiten
 subplot(ceil(length(h_vec)/2)+1,2,...
     (length(h_vec)+2):(2*ceil(length(h_vec)/2)+2))
 plot(diff(t));
 axis([0 length(t)-1 min(diff(t)) max(diff(t))])
 title('Schrittweiten in ode23') 
 %Speichern des Diagramms
 print('-depsc','pl4_3.eps')
end

%Ein Schritt des Verfahrens von Heun
function[y_neu]=heun_step(fct,t,y,h)
% !!!!! Implementiere hier einen Schritt des Verfahrens von Heun  !!!!!
end

function[f_]=yprime(t,y)
.....
end
