% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% --- Programmieraufgabe 1 ---
% Namen der Gruppenmitglieder:
% 1. ...
% 2. ...
% 3. ...
% 
% Alle Dateien muessen bis zum 31.10.2019 hochgeladen werden (siehe
% ISIS-Webseite)
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
close all;
N_samples=[10,50,100,500]; % Test multiple nbs of steps

%% 2.a)
% Set the problem 
t0  = 0;
y0  = 1;
a   = 1.5;
f   = @(t,y) -tan(t).*y;      % Right hand side of the ODE y' = f(t,y)
y = @(t) cos(t);

% Preparation of the plots
figure(1)
title('dy/dt=-tan(t)y');
xlabel('t');
ylabel('y');
hold on;
legendLabels=cell(2*length(N_samples)+1, 1); % Initialization
err_EE_a=zeros(length(N_samples),1);
err_Heun_a=err_EE_a;
h_a = zeros(length(N_samples), 1);

for i=1:length(N_samples)
    N = N_samples(i);
    [h,t,u] = expeuler_loesung(f, t0,y0, N, a);
    err_EE_a(i)=abs(u(end)-y(t0+a));
    h_a(i)=h;
    
    figure(1) % Plot the current solution y(t)
    plot(t, u, '-.o', 'MarkerSize', 3);
    legendLabels{i}=['N=',num2str(N),' EE'];
end
for i=1:length(N_samples)
    N = N_samples(i);
    [h,t,u] = heun_loesung(f, t0,y0, N, a);
    err_Heun_a(i)=abs(u(end)-y(t0+a));
    
    figure(1) % Plot the current solution y(t)
    plot(t, u, '-.o', 'MarkerSize', 3);
    legendLabels{length(N_samples)+i}=['N=',num2str(N),', Heun'];
end
figure(1)
plot(t, y(t), '-g')
legendLabels{end}=['Analytische Loesung'];
legend(legendLabels, 'Location', 'NorthWest');

%% 2.b)
% Set the problem 
t0  = 0;
y0  = [sin(1); cos(1)];
a   = 2;
f   = @(t,y) [0, 1; -exp(2*t), 1]*y;      % Right hand side of the ODE y' = f(t,y)
y = @(t) [sin(exp(t)); cos(exp(t)).*exp(t)];

% Preparation of the plots
figure(21)
title('2.b) erster Eintrag');
xlabel('t');
ylabel('y_1');
hold on;
figure(22)
title('2.b) zweiter Eintrag');
xlabel('t');
ylabel('y_2');
hold on;
legendLabels=cell(2*length(N_samples)+1, 1); % Initialization

err_EE_b=zeros(length(N_samples),1);
err_Heun_b=err_EE_b;
h_b = zeros(length(N_samples), 1);

for i=1:length(N_samples)
    N = N_samples(i);
    [h,t,u] = expeuler_loesung(f, t0, y0, N, a);
    err_EE_b(i)=norm(u(:,end)-y(t0+a),2);
    h_b(i)=h;
    
    figure(21) % Plot the current solution y(t)
    plot(t, u(1,:), '-.o', 'MarkerSize', 3);
    legendLabels{i}=['N=',num2str(N),' EE'];
    
    figure(22) % Plot the current solution y(t)
    plot(t, u(2,:), '-.o', 'MarkerSize', 3);
end
for i=1:length(N_samples)
    N = N_samples(i);
    [h,t,u] = heun_loesung(f, t0, y0, N, a);
    err_Heun_b(i)=norm(u(:,end)-y(t0+a),2);
    
    figure(21) % Plot the current solution y(t)
    plot(t, u(1,:), '--x', 'MarkerSize', 3);
    legendLabels{length(N_samples)+i}=['N=',num2str(N),', Heun'];
    
    figure(22) % Plot the current solution y(t)
    plot(t, u(2,:), '--x', 'MarkerSize', 3);
end
yt=y(t);
figure(21)
plot(t, yt(1,:), '-g')
legendLabels{end}='Analytische Loesung';
legend(legendLabels, 'Location', 'NorthWest');

figure(22)
plot(t, yt(2,:), '-g')
legendLabels{end}='Analytische Loesung';
legend(legendLabels, 'Location', 'NorthWest');

%% 2.c)
figure(3)
loglog(h_a, err_EE_a, 'x');
hold on;
loglog(h_a, err_Heun_a, 'o');
title('Konvergenzdiagramm fuer die Aufgabe a)')
xlabel('N')
ylabel('err')
legend('EE', 'Heun', 'Location', 'NorthWest');

q_EE_a=polyfit(log(h_a), log(err_EE_a),1);
disp(['a) Slope of the regression line (EE): ', num2str(q_EE_a(1))]);

q_Heun_a=polyfit(log(h_a), log(err_Heun_a),1);
disp(['a) Slope of the regression line (Heun): ', num2str(q_Heun_a(1))]);

figure(4)
loglog(h_b, err_EE_b, 'x');
hold on;
loglog(h_b, err_Heun_b, 'o');
title('Konvergenzdiagramm fuer die Aufgabe b)')
xlabel('h')
ylabel('err')
legend('EE', 'Heun');

q_EE_b=polyfit(log(h_b), log(err_EE_b),1);
disp(['b) Slope of the regression line (EE): ', num2str(q_EE_b(1))]);

q_Heun_b=polyfit(log(h_b), log(err_Heun_b),1);
disp(['b) Slope of the regression line (Heun): ', num2str(q_Heun_b(1))]);

% Terminate the sentence below.
disp(['=>The convergence order of EE is ... ']);
disp(['=>The convergence order of Heuns method is ... ']);






