% m27_RK4convergence.m - Solve the scalar IVP  u' = sin(40tu), u(0)=1

% First for comparison we solve it by ode45:
  FS = 'fontsize'; IN = 'interpreter'; LT = 'latex'; CO =  'color';
  tmax = 1; u0 = 1;
  opts = odeset('reltol',1e-8,'abstol',1e-8);
  f = @(t,u) sin(40*t*u);
  [t,u] = ode45(f,[0 tmax],1,opts);
  hold off, plot(t,u,'-r','linewidth',2)
  ylim([.7 1.2]), grid on, hold on
  title('Exact solution',FS,18,CO,'r')
  exact = u(end);

% Now we solve it by 4th-order RK with various fixed step sizes k:
  pause
  err = [];
  for k = .5.^(1:7)
    nsteps = round(tmax/k);
    u = u0; uu = u; tt = 0;
    for n = 0:nsteps-1
      t = n*k;
      a = k*f(t,u);
      b = k*f(t+k/2,u+a/2);
      c = k*f(t+k/2,u+b/2);
      d = k*f(t+k,u+c);
      u = u + (a+2*b+2*c+d)/6;          % <-- Runge-Kutta formula
      uu = [uu; u]; tt = [tt; t+k];
    end
    if k==.5
       p = plot(tt,uu,'.-k','markersize',18,'linewidth',.6); grid on
    else
       set(p,'xdata',tt,'ydata',uu)
    end
    err = [err; abs(uu(end)-exact)];
    title(['$k = ' num2str(k) '$'],FS,18,CO,'k',IN,LT);
    pause
  end

% Finally we plot the errors as a function of step size k:
  hold off, semilogy(1:7,abs(err),'.-b','markersize',18)
  xlabel('| log_2 k |',FS,16)
  ylabel('error at $t=1$',FS,16,IN,LT), grid on, hold on
  semilogy(3:7,10*2.^(-4*(3:7)),'--r','linewidth',2) 
  text(4.3,6e-6,'$k^4$',FS,20,CO,'r',IN,LT)
  title('Fourth-order convergence',FS,20)