%====================================================================
% Initial value problem for a single 1st-order ODE
% Call Methods: Euler, Predictor-Corrector, Runge-Kutta 4th order
%====================================================================
clear;
%open a file if needed
% fileID = fopen('results.txt','w');
% initial conditions
xi = 1.00;
ti = 0.00;
% "time" step and max time
dt =   0.01;
tmax = 10.0;
% select method: key=1 Euler, key=2 predictor-corrector, key=3 RK 4th order
key = 3;
%select function
dx = @f;

%preparing to plot
i    =  1;
t (i)= ti;
x1(i)= xi;

% print the header and initial conditions
fprintf(' Solving 1st order single ODE \n')
switch key
    case 1
        fprintf(' Simple Euler method \n')
    case 2
        fprintf(' Predictor-Corrector \n')
    case 3
        fprintf(' Runge-Kutta 4th order \n')
end
fprintf('     t         x(t)   \n')
formatSpec1 = '%10.5f %10.5f \n';
%fprintf(fileID,formatSpec1,ti,xi)
fprintf(formatSpec1,ti,xi)

while(ti < tmax)
  tf = ti + dt;
  i = i+1;
  if(key==1) 
      [xf] = euler1(dx,ti,xi,tf);
  end
  if(key==2) 
      [xf] = euler2(dx,ti,xi,tf);
  end
  if(key==3) 
      [xf] = RK4D11(dx,ti,xi,tf);
  end
%  xtrue = exp(-1.0*tf);
  fprintf(formatSpec1,tf,xf); 
%  fprintf(fileID,'%10.5f %10.5f %10.5f\n',tf,xf,xtrue);
% to plot
  t(i)  = tf;
  x1(i) = xf;
  
  ti = tf;
  xi = xf;
end %while
plot(t,x1,'r')


function f=f(t,x)
%========
% function for dx/dt = f(t,x)
%========
%    f = (-1.0)*x+t*0.0;  
   f = (-1.0)*x*cos(t)+0.5*sin(2*x);
% f = t*x*x;
end

function [xf] = euler1(dx,ti,xi,tf)
%==================================================
% Solution of a single 1st order ODE dx/dt=f(x,t)
% Method:  Simple Euler
% Alex G. October 2020
%--------------------------------------------------
% input ...
% dx(t,x)- function dx/dt (supplied by a user)
% ti - initial time
% xi  - initial position
% tf  - time for a solution
% output ...
% xf  - solution at point tf
%==================================================
xf = xi + dx(ti,xi)*(tf-ti);
end

function [xf] = euler2(dx,ti,xi,tf)
%==================================================
% Solution of a single 1st order ODE dx/dt=f(x,t)
% Method:  Euler: predictor-corrector
% Alex G. October 2020
%--------------------------------------------------
% input ...
% dx(t,x)- function dx/dt (supplied by a user)
% ti - initial time
% xi  - initial position
% tf  - time for a solution
% output ...
% xf  - solution at point tf
%==================================================
xf = xi + dx(ti,xi)*(tf-ti);
xf = xi + (dx(ti,xi)+dx(tf,xf))*0.5*(tf-ti);
end

function [xf] = RK4D11(dx,ti,xi,tf)
%==================================================
% Solution of a single 1st order ODE dx/dt=f(x,t)
% Method:  Runge-Kutta 4th order
% Alex G. October 2020
%--------------------------------------------------
% input ...
% dx(t,x)- function dx/dt (supplied by a user)
% ti - initial time
% xi  - initial position
% tf  - time for a solution
% output ...
% xf  - solution at point tf
%==================================================
h  = tf-ti;

k1 = h*dx(ti,xi);
k2 = h*dx(ti+h/2.0,xi+k1/2.0);
k3 = h*dx(ti+h/2.0,xi+k2/2.0);
k4 = h*dx(ti+h,xi+k3);

xf = xi + (k1 + 2.0*(k2+k3) + k4)/6.0;
end