%{
Solving 1st-order ODE
Method: RKF45 with step-size control
Alex G. March 2022
%}

clear;
clc;
close;

dx = @f;
fa = @(t) atan(100*t);     % analytic solution for function f (if availble)
KeyA = 0;    % 1 - if comparing with analytic solution, 0 - if not

tolerance = 1.e-6;
StepsMax  = 1000;

% 1: initial condition
% ti = 0.00;
% xi = 1.00;
ti = -1.0;
xi = atan(-100);

% 2: "time" step and max time
dt =   0.5;
tmax = 1.0;

% 3: prepare to plot an analytic solution (if available)
if KeyA == 1
  a = ti;
  b = tmax+dt;
  Nplot = 100;
  h = (b-a)/Nplot;
  for k = 1:Nplot+1
      xa(k) = a+ (k-1)*h;
      ya(k) = fa(xa(k));
  end
end

% 4: getting ready for the main loop
i    = 1;
step = 1;
t(i) = ti;
x(i) = xi;

% 5: the MAIN part
while (ti < tmax)
  tf = ti + dt;
  [xf,er5] = RKF51(dx,ti,xi,tf);
  er5 = abs(er5);
  if (er5 <= tolerance)
    i = i + 1;
    x(i) = xf;
    t(i) = tf;
    ti = tf;
    xi = xf;
    dt = min(0.9*dt*(abs(tolerance/er5)).^(1/5),2.0*dt); 
  else
    dt = max(0.9*dt*(abs(tolerance/er5)).^(1/5),0.2*dt);
  end
  step = step + 1;
  if step > StepsMax
    fprintf(' Max number of steps exceeded at tf = %8.4f \n',tf)
    break
  end
end

% figure (1)
figure('DefaultAxesFontSize',14)
if KeyA == 0
    plot(t,x,'b-o','LineWidth',1.0)
    legend('RKF45')
  elseif KeyA == 1
    plot(t,x,'bo',xa,ya,'r','LineWidth',1.0)
    legend('RKF45','analytic solution')
end
xlabel('t')
ylabel('x(t)')

grid

function f=f(t,x)
%========
% function for dx/dt = f(t,x)
%========
%f = 0.0-x+0*t;
%f = x+exp(-t);
%f = 0.5*sin(2*x) - x*cos(t);
f = 100/(1+10000.*t.^2);
end

function [xf,er5] = RKF51(dx,ti,xi,tf)
%==================================================
% Solution of a single 1st order ODE dx/dt=f(x,t)
% Method:  Runge-Kutta-Fehlberg RKF45
% Alex G. March 2022
%--------------------------------------------------
% 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
% er511 - error estimate
%==================================================
h  = tf-ti;
k1 = h*dx(ti,xi);
k2 = h*dx(ti+h/4.0,  xi+k1/4.0);
k3 = h*dx(ti+3*h/8,  xi+3*k1/32+9*k2/32);
k4 = h*dx(ti+12*h/13,xi+1932*k1/2197-7200*k2/2197+7296*k3/2197);
k5 = h*dx(ti+h,      xi+439*k1/216-8*k2+3680*k3/513-845*k4/4104);
k6 = h*dx(ti+h/2,    xi-8*k1/27+2*k2-3544*k3/2565+1859*k4/4104-11*k5/40);
xf = xi + 16*k1/135+6656*k3/12825+28561*k4/56430-9*k5/50+2*k6/55;
er5 = k1/360-128*k3/4275-2197*k4/75240+k5/50+2*k6/55;
end