%% Numerical implemations - Section 2.4, 2.5, 2.6
% One program to do it all: Euler, Heun and Runge Kutta
% Interactive all inputs are given on the command window - nothing to
% change in the main code
% Use:  numerics; - run the code, doesn't save data
%       [xE,yE,~,~,~,~]=numerics; - run the code, save the Euler data
%       ....
%       [xE,yE,xH,yH,xRK,yRK]=numerics; - run the code, save all data

% In subsequent releases, inline is going to be removed. Then replaced by
% syms f(x,y) before the input call and delete 's' in input.

function [xE,yE,xH,yH,xRK,yRK]=numerics

% Initialization - Interactive via command window
method = input('Enter the method you want to use: Euler, Heun, RK, all: ','s');
fprintf('\n')
output = input('Enter the type of output: table, plot, both: ','s');
fprintf('\n')
f = input('Enter the right-hand side of the ODE as a function of y and x: ','s');
f = inline(f,'x','y');
x0 = input('Enter the initial x value x0: ');
y0 = input('Enter the initial y value y0: ');
fprintf('\n')
h = input('Enter the step size h: ');
querry = input('Do you know the final x value? [Y/N]: ','s');
if querry=='Y'
    xf = input('Enter the final x value: ');
    N = (xf-x0)/h;
else
    N = input('Enter the number of steps: ');
    xf = x0+N*h;
end
fprintf('\n')
exact = input('Do you know the exact solution? [Y/N]: ','s');
if exact=='Y'
    yexact = input('Enter the exact solution: ','s');
    yexact = inline(yexact,'x');
    err = input('Do you want to calculate the global error at the final point? [Y/N]: ','s');
end
if exact == 'N'
    err = 'N';
end
fprintf('\n')

% Main code: run one or all methods and output results in a table or plot
if strcmp(method,'Euler') || strcmp(method,'all')
    [xE,yE] = Euler(f,x0,y0,h,N);
    if strcmp(method,'Euler')
        xH = [];
        yH = [];
        xRK = [];
        yRK = [];
    end
    if strcmp(output,'plot') || strcmp(output,'both')
        plot(xE,yE,'LineWidth',1.5)
        hold on
    end
    if strcmp(output,'table') || strcmp(output,'both') 
        fprintf('Euler:\n')
        X = ['x ' num2str(xE,'%10.6f')];
        Y = ['y ' num2str(yE,'%10.6f')];
        disp(X)
        fprintf('\n')
        disp(Y)
        fprintf('\n')
    end
end
if strcmp(method,'Heun') || strcmp(method,'all')
    [xH,yH] = Heun(f,x0,y0,h,N);
    if strcmp(method,'Heun')
        xE = [];
        yE = [];
        xRK = [];
        yRK = [];
    end
    if strcmp(output,'plot') || strcmp(output,'both')
        plot(xH,yH,'r','LineWidth',1.5)
    end
    if strcmp(output,'table') || strcmp(output,'both') 
        fprintf('Heun:\n')
        X = ['x ' num2str(xH,'%10.6f')];
        Y = ['y ' num2str(yH,'%10.6f')];
        disp(X)
        fprintf('\n')
        disp(Y)
        fprintf('\n')
    end
end
if strcmp(method, 'RK') || strcmp(method,'all')
    [xRK,yRK] = RK(f,x0,y0,h,N);
    if strcmp(method,'RK')
        xH = [];
        yH = [];
        xE = [];
        yE = [];
    end
    if strcmp(output,'plot') || strcmp(output,'both')
        plot(xRK,yRK,'m','LineWidth',1.5)
    end
    if strcmp(output,'table') || strcmp(output,'both')
        fprintf('Runge Kutta:\n')
        X = ['x ' num2str(xRK,'%10.6f')];
        Y = ['y ' num2str(yRK,'%10.6f')];
        disp(X)
        fprintf('\n')
        disp(Y)
        fprintf('\n')
    end                
end
if (strcmp(output,'plot')||strcmp(output,'both')) && strcmp(exact,'Y')
    plot(x0:h:xf,yexact(x0:h:xf),'k','Linewidth',1.5)
end
if (strcmp(output,'table')||strcmp(output,'both')) && strcmp(exact,'Y')
    fprintf('Exact solution:\n')
    xp = x0+(0:N)*h;
    X = ['x ' num2str(xp,'%10.6f')];
    Y = ['y ' num2str(yexact(xp),'%10.6f')];
    disp(X)
    fprintf('\n')
    disp(Y)
    fprintf('\n')
end

if strcmp(err,'Y')
    if (strcmp(method,'Euler')||strcmp(method,'all'))
        errorE = abs(yexact(xf)-yE(end));    
        relerrorE = errorE/abs(yexact(xf))*100;
        fprintf('Euler global error %10.6f\n',errorE);
        fprintf('Euler relative error %10.6f%%\n', relerrorE);
    end
    if (strcmp(method,'Heun')||strcmp(method,'all'))
        errorH = abs(yexact(xf)-yH(end));
        relerrorH = errorH/abs(yexact(xf))*100;
        fprintf('Heun global error %10.6f\n',errorH);
        fprintf('Heun relative error in %10.6f%%\n', relerrorH);
    end
    if (strcmp(method,'RK')||strcmp(method,'all'))
        errorRK = abs(yexact(xf)-yRK(end));
        relerrorRK = errorRK/abs(yexact(xf))*100;
        fprintf('Runge Kutta global error %10.6f\n',errorRK);
        fprintf('Runge Kutta relative error in %10.6f%%\n', relerrorRK);
    end
end

if (strcmp(output,'plot')||strcmp(output,'both'))
    xlabel('x','FontSize',12)
    ylabel('y','FontSize',12)
    set(gcf,'color','w')
    set(gca,'FontSize',12)
    set(get(gca,'YLabel'),'Rotation',0)
    if strcmp(method,'Euler')
        legend('Euler','Location','NorthWest')
    elseif strcmp(method,'Heun')
        legend('Heun','Location','NorthWest')
    elseif strcmp(method,'RK')
        legend('RK','Location','NorthWest')
    else
    legend('Euler','Heun','RK','Location','NorthWest')
    end
hold off
end


%% Euler method
% Solve the IVP y'(x)=f(x,y) with y(y0)=x0. The function f is given as a
% symbolic function.
% Given: step size h and number of steps. 

function [xE,yE] = Euler(f,x0,y0,h,N)

xE = zeros(1,N+1); yE = zeros(1,N+1);   % Initialization
xE(1) = x0; yE(1) = y0;                 % Initial condition

for i=1:N                               % For loop over the steps
    xE(i+1) = xE(i)+h;                  % Next x
    yE(i+1) = yE(i)+h*f(xE(i),yE(i));   % Euler iteration - next y
end


%% Heun method
% Solve the IVP y'(x)=f(x,y) with y(y0)=x0. The function f is given as a
% symbolic function.
% Given: step size h and number of steps. 

function [xH,yH] = Heun(f,x0,y0,h,N)

xH = zeros(1,N+1); yH = zeros(1,N+1);   % Initialization 
xH(1) = x0; yH(1) = y0;                 % Initial condition

for i=1:N                               % For loop over the steps
    k1 = f(xH(i),yH(i));                % k1
    k2 = f(xH(i)+h,yH(i)+h*k1);         % k2
    k = (k1+k2)/2;                      % k
    xH(i+1) = xH(i)+h;                  % Next x
    yH(i+1) = yH(i)+h*k;                % Heun iteration - next y
end


%% RK method
% Solve the IVP y'(x)=f(x,y) with y(y0)=x0. The function f is given as a
% symbolic function.
% Given: step size h and number of steps. 

function [xRK,yRK] = RK(f,x0,y0,h,N)

xRK = zeros(1,N+1); yRK = zeros(1,N+1); % Initialization
xRK(1) = x0; yRK(1) = y0;               % Initial condition

for i=1:N                               % For loop over the steps
    k1 = f(xRK(i),yRK(i));              % k1
    k2 = f(xRK(i)+h/2,yRK(i)+h*k1/2);   % k2       
    k3 = f(xRK(i)+h/2,yRK(i)+h*k2/2);   % k3
    k4 = f(xRK(i)+h,yRK(i)+h*k3);       % k4
    k = (k1+2*k2+2*k3+k4)/6;            % k
    xRK(i+1) = xRK(i)+h;                % Next x
    yRK(i+1) = yRK(i)+h*k;              % Runge Kutta iteration - next y
end
