%% MCHI, cv4, file 1 (cv4_1) - Numerical solution to Cauchy problem
%
% Author:       Martin Isoz
% Organisation: ICT Prague
% Date:         09. 10. 2014
%
% License: This code is published under MIT License, please do not abuse
% it.
%

% Note: dxdt = f(x) ... autonomous system

%% Caller 'script'
function caller
close all

% solving problem #1
Y0      = 1;
TSPAN   = [0 1];
H       = 0.1;

tstart = tic;
[TOUT,YOUTE] = eulerM(@problem1,TSPAN,Y0,H);
timeE = toc(tstart);
tstart = tic;
[~,YOUTRK4]  =   rk4M(@problem1,TSPAN,Y0,H);
timeRK4 = toc(tstart);

figure
plot(TOUT,exp(TOUT),'k-');hold on;plot(TOUT,YOUTE,'ko');plot(TOUT,YOUTRK4,'ks');
axis tight
title('Solution of problem 1','FontSize',18,'FontWeight','Bold')
legend('Exact',['Euler t = ' num2str(timeE) ' s'],['RK4 t = ' num2str(timeRK4) ' s'],'Location','Best');

% solving problem #2
Y0      = [2/3 1];
TSPAN   = [0 1];
H       = 0.1;

[TOUT,YOUTE] = eulerM(@problem2,TSPAN,Y0,H);
[~,YOUTRK4]  =   rk4M(@problem2,TSPAN,Y0,H);

figure
plot(TOUT,exactSol2(TOUT));hold on;plot(TOUT,YOUTE,'o');plot(TOUT,YOUTRK4,'s');
legend('xExact','yExact','xEuler','yEuler','xRK4','yRK4','Location','Best')
% plot(TOUT,exp(TOUT),'k-');hold on;plot(TOUT,YOUTE,'ko');plot(TOUT,YOUTRK4,'ks');
axis tight
title('Solution of problem 2','FontSize',18,'FontWeight','Bold')

% solving problem #3
% this problem will be solved only using RK4M solver but for the different
% set of initial conditions

% create a matrix for different initial conditions
Y0M = [0.2:0.1:1;0.2:0.1:1]';                                                 %create matrix 2xN than transpose it to Nx2

% prepare the rest of inputs for solver
TSPAN   = [0 10];
H       = 0.1;

% prepare parameters for the model
alpha = 1;beta = 1;gamma = 1;delta = 1;
pars  = [alpha beta gamma delta];                                           %same structure as in problem3 function

% prepare figure window
figure;hold on

% solve the problem
for i = 1:numel(Y0M)/2
    Y0 = Y0M(i,:);                                                          %I will start in point (0.2,0.2),(0.4,0.4)...
    [~,YOUT]  =   rk4M(@(t,x) problem3(t,x,pars),TSPAN,Y0,H);               %this time I am not interested in dependent variable
    plot(YOUT(:,1),YOUT(:,2),'.','Color',rand(1,3));                        %plot trajectory in the phase space
end

% format the graphics
title('Trajectories of problem3 in plane x_1-x_2','FontSize',18,'FontWeight','Bold')
xlabel('x_1');ylabel('x_2');

end

%% Methods
function [TOUT,YOUT] = eulerM(ODEFUN,TSPAN,Y0,H)                            %with tspan we input time step (not adaptive)

% calculate needed # of steps
nSteps = norm(diff(TSPAN))./H+1;                                            %I will not end exactly in TEND

% allocate variables
TOUT = zeros(nSteps,1); YOUT = zeros(nSteps,numel(Y0));

% save initial condition
TOUT(1) = min(TSPAN);YOUT(1,:) = reshape(Y0,1,numel(Y0));

for i = 2:nSteps
    TOUT(i) = TOUT(i-1) + H;
    YOUT(i,:) = YOUT(i-1,:) + H*reshape(feval(ODEFUN,TOUT(i-1),YOUT(i-1,:)),1,numel(Y0));
end

end


function [TOUT,YOUT] = rk4M(ODEFUN,TSPAN,Y0,H)                              %with tspan we input time step (not adaptive)
% calculate needed # of steps
nSteps = norm(diff(TSPAN))./H+1;                                            %I will not end exactly in TEND

% allocate variables
TOUT = zeros(nSteps,1); YOUT = zeros(nSteps,numel(Y0));

% save initial condition
TOUT(1) = min(TSPAN);YOUT(1,:) = reshape(Y0,1,numel(Y0));

for i = 2:nSteps
   TOUT(i) = TOUT(i-1) + H; 
   k1 = H*reshape(feval(ODEFUN,TOUT(i-1),YOUT(i-1,:)),1,numel(Y0));
   k2 = H*reshape(feval(ODEFUN,TOUT(i-1)+0.5*H,YOUT(i-1,:)+0.5*k1),1,numel(Y0));
   k3 = H*reshape(feval(ODEFUN,TOUT(i-1)+0.5*H,YOUT(i-1,:)+0.5*k2),1,numel(Y0));
   k4 = H*reshape(feval(ODEFUN,TOUT(i-1)+1*H,YOUT(i-1,:)+1*k3),1,numel(Y0));
   YOUT(i,:) = YOUT(i-1,:) + 1/6*(k1+k4) + 1/3*(k2+k3);
end
end

%% Problem definitions
% first equation to be solved (x in R)
function dxdt = problem1(~,x)                                               %I do not need to input t -> ~
dxdt = x;
end

% second equation to be solved (x in R2)
function dxdt = problem2(~,x)                                               %%I do not need to input t -> ~

dxdt = zeros(numel(x),1);                                                   %allocate space for output (general)

dxdt(1) = x(1) + x(2).^2;                                                   %problem specific
dxdt(2) = -x(2);                                                            %problem specific
end

% exact solution to the problem 2
function xExact = exactSol2(t)

xExact = zeros(numel(t),2);
xExact(:,1) = exp(t) - exp(-2*t)./3;
xExact(:,2) = exp(-t);
end

function dxdt = problem3(~,x,pars)
% Note: This is a SODE defining the classical predator-prey model
% (Lotka-Voltera equations)

dxdt = zeros(numel(x),1);

% unwrap parameters
alpha = pars(1);
beta  = pars(2);
gamma = pars(3);
delta = pars(4);

% calculate the residues
dxdt(1) = x(1)*(alpha - beta*x(2));
dxdt(2) = -x(2).*(gamma - delta*x(1));
end