function aufg18

% initial value in polar coordinates
r_0     = 2;
theta_0 = 0;
u_ini   = r_0*[cos(theta_0); sin(theta_0)]; % in Euclidean coordinates

% time grid
t_0   = 0;
t_end = 10000;
% h = 0.1; 
h = 0.2;

t_grid = t_0 : h : t_end;

% error bounds
c_eul = 1;
c_rk  = 0.1;

% animation? 0 == no, 1 == yes
animation = 1;

% calculating the analytic solution
u_sol = u_ana(t_grid,r_0,theta_0);

% calculating the numerical solutions
f = @(t,u) [-u(2) + u(1)*(1 - u(1)^2 - u(2)^2);...
    u(1) + u(2).*(1-u(1)^2 - u(2)^2)] + 0*t;
u_eul = euler(t_grid,u_ini,f);
u_rk  = rkclassic(t_grid,u_ini,f);

err_euler = 0;
err_rk    = 0;
t_err_eul = 0;
t_err_rk  = 0;

% computing the error and error time for the Euler scheme
for j = 1:length(t_grid)
    err_euler = norm(u_eul(:,j) - u_sol(:,j));
    t_err_eul = h*(j-1);
    if err_euler >= c_eul
        break
    end    
end

% computing the error and error time for the RK scheme
for j = 1:length(t_grid)    
    err_rk    = norm(u_rk(:,j)  - u_sol(:,j));
    t_err_rk = h*(j-1);
    if err_rk >= c_rk
        break
    end
end

fprintf('\nAufgabe 18\n')
fprintf(['stepsize h = ' num2str(h) '\n'])
fprintf(['Euler method: error = ' num2str(err_euler) ' at time = '...
    num2str(t_err_eul) '\n'])
fprintf(['  classic RK: error = ' num2str(err_rk) ' at time = '...
    num2str(t_err_rk) '\n\n']) 

% generating the animation
if animation == 1  

figure(2)
axis([-2 2 -2 2])
j = 1;
err_euler = 0;
err_rk = 0;
while j <= length(t_grid)
    k  = round(5/h);
    j0 = max(1,j-k);
    n  = 0.01;
    
    % plotting the analytic solution
    plot(u_sol(1,j0:j),u_sol(2,j0:j),'r')
    axis([-2 2 -2 2])
    hold on
    plot(u_sol(1,j),u_sol(2,j),'r*')
    title(['t = ' num2str((j-1)*h)])
    xlabel('u_1')
    ylabel('u_2')
    
    % plotting the Euler approximation
    if err_euler < c_eul
        plot(u_eul(1,j0:j),u_eul(2,j0:j),'k')
        plot(u_eul(1,j),u_eul(2,j),'k*')
        err_euler = norm(u_eul(:,j) - u_sol(:,j));
        if err_euler >= c_eul
            n=0.01;
            pause
            %j = j + round(0.3*t_end/h);
        end
    end
    
    % plotting the RK-approximation
    if err_rk < c_rk
        plot( u_rk(1,j0:j), u_rk(2,j0:j),'b')
        plot( u_rk(1,j),u_rk(2,j),'b*')
        err_rk = norm(u_rk(:,j) - u_sol(:,j));
        t_err_rk = h*(j-1);
        title(['t = ' num2str(t_err_rk)])
        xlabel('u_1')
        ylabel('u_2')
        if err_rk >= c_rk
            break
        end
    end
    hold off
    
    pause(n)
    j = j+3;
end  

end

end

% auxiliary functions

% analytic solution
function out = u_ana(t,r_0,theta_0)
% t = time grid
% r_0, theta_0 = polar coordinates of the initial value
r = sqrt(1 ./ ( (1-r_0^2)*exp(-2*t)/r_0^2 + 1) );
theta = t + theta_0;
out = [r.*cos(theta);r.*sin(theta)];

end

% Euler method
function y = euler(t,u_ini,f)
m = length(u_ini); % dimension of the system
l = length(t);     % number of points in the grid

y = zeros(m,l);    % initialising the variable for the results
y(:,1) = u_ini;    % setting the initial value
% Euler scheme
for i = 1:(l-1)
    h = t(i+1) - t(i);
    y(:,i+1) = y(:,i) + h*f(t(i),y(:,i));
end
end

% classic Runge-Kutta method
function y = rkclassic(t,u_ini,f)
%  input: t     ~ vector of time points
%         u_ini ~ initial value at time t(1)
%         f     ~ function handle for the ODE u'(t) = f(t,u(t))
% output: y     ~ numerical approximation of the ODE, y(n) ~ u(t(n))
    
m = length(u_ini); % dimension of the system
l = length(t);     % number of points in the grid

% RK classic coefficients
alpha = [0 1/2 1/2 1];
beta  = [0 0 0 0; 1/2 0 0 0; 0 1/2 0 0; 0 0 1 0];
gamma = [1 2 2 1]/6;

y = zeros(m,l);    % initialising the variable for the results
y(:,1) = u_ini;    % setting the initial value
% Heun scheme
for i = 1:(l-1)
    h = t(i+1) - t(i);
    k = zeros(m,length(alpha)); % memory for the k_j coefficient functions
    % calculating the k_j, compare with (2.3)
    for j = 1:length(alpha)
        k(:,j) = f(t(i) + alpha(j)*h,y(:,i)+h*k*beta(j,:)');
    end
    y(:,i+1) = y(:,i) + h*k*gamma';
end
end
