clc
clf
clear all

%    RK4 method

% Solve dy/dt = f(t,y). In general it can be a function of both 
% variables t and y. If your function is only a function of t then
% you will need to add a 0*y to your function.

%  define f(t,y)
   fcnstr='y-t^2+1' ;
   f=inline(fcnstr,'t','y') ;

% t0, initial time

   t0=0 ;

% y0, corresponding value of y at t0

   y0=0.5;

% tf, upper boundary of time t (end time). 

   tf=2;

% n, number of steps to take

   n=5;

%**********************************************************************

% Displays title information
disp(sprintf('\n\nThe 4th Order Runge-Kutta Method of Solving Ordinary Differential Equations'))
h=(tf-t0)/n ;
disp(sprintf('     h = ( tf - t0 ) / n '))
disp(sprintf('       = ( %g - %g ) / %g ',tf,t0,n))
disp(sprintf('       = %g',h))

ta(1)=t0 ;
ya(1)=y0 ;


for i=1:n
  disp(sprintf('\nStep %g',i))
  disp(sprintf('-----------------------------------------------------------------'))

% Adding Step Size
  ta(i+1)=ta(i)+h ;

% Calculating k1, k2, k3, and k4
  k1 = f(ta(i),ya(i)) ;
  k2 = f(ta(i)+0.5*h,ya(i)+0.5*k1*h) ;
  k3 = f(ta(i)+0.5*h,ya(i)+0.5*k2*h) ;
  k4 = f(ta(i)+h,ya(i)+k3*h) ;
     
% Using 4th Order Runge-Kutta formula
  ya(i+1)=ya(i)+1/6*(k1+2*k2+2*k3+k4)*h ;

  disp('1) Find k1 and k2 using the previous step information.')
  disp(sprintf('     k1 = f( x%g , y%g )',i-1,i-1))
  disp(sprintf('        = f( %g , %g )',ta(i),ya(i)))
  disp(sprintf('        = %g\n',k1))

  disp(sprintf('     k2 = f( x%g + 0.5 * h , y%g + 0.5 * k1 * h )',i-1,i-1))
  disp(sprintf('        = f( %g + 0.5 * %g , %g + 0.5 * %g * %g)',ta(i),h,ya(i),k1,h))
  disp(sprintf('        = f( %g , %g )',ta(i)+0.5*h,ya(i)+0.5*k1*h))
  disp(sprintf('        = %g\n',k2))

  disp(sprintf('     k3 = f( x%g + 0.5 * h , y%g + 0.5 * k2 * h )',i-1,i-1))
  disp(sprintf('        = f( %g + 0.5 * %g , %g + 0.5 * %g * %g)',ta(i),h,ya(i),k2,h))
  disp(sprintf('        = f( %g , %g )',ta(i)+0.5*h,ya(i)+0.5*k2*h))
  disp(sprintf('        = %g\n',k3))

  disp(sprintf('     k4 = f( x%g + h, y%g + k3*h)',i-1,i-1))
  disp(sprintf('        = f( %g , %g )',ta(i)+h,ya(i)+k3*h))
  disp(sprintf('        = %g\n',k1))

  disp(sprintf('2) Apply the Runge-Kutta 4th Order method to estimate y%g',i))
  disp(sprintf('     y%g = y%g + 1/6*( k1 + 2*k2 + 2*k3 +k4) * h',i,i-1))
  disp(sprintf('        = %g + %g * %g * %g',ya(i),1/6,(k1+2*k2+2*k3+k4),h))
  disp(sprintf('        = %g\n',ya(i+1)))
  disp(sprintf('   at t%g = %g',i,ta(i+1)))
end


disp(sprintf('\n\n********************************Results**********************************'))

% The following finds what is called the 'Exact' solution
tspan = [t0 tf];
[t,y]=ode45(f,tspan,y0);
[yfi dummy]=size(y);
yf=y(yfi);

% Plotting the Exact and Approximate solution of the ODE.
hold on
xlabel('x');ylabel('y');
title('Exact and Approximate Solution of the ODE by the 4th Order Runge-Kutta Method');
plot(t,y,'--','LineWidth',2,'Color',[0 0 1]);            
plot(ta,ya,'-','LineWidth',2,'Color',[0 1 0]);
legend('Exact','Approximation');

disp('The figure window that now appears shows the approximate solution as ') 
disp('piecewise continuous straight lines. The blue line represents the exact')
disp('solution. In this case ''exact'' refers to the solution obtained by the')
disp(sprintf('Matlab function ode45.\n'))

disp('Note the approximate value obtained as well as the true value and ')
disp('relative true error at our desired point x = xf.')

disp(sprintf('\n   Approximate = %g',ya(n+1))) 
disp(sprintf('   Exact       = %g',yf))
disp(sprintf('\n   True Error = Exact - Approximate')) 
disp(sprintf('              = %g - %g',yf,ya(n+1)))
disp(sprintf('              = %g',yf-ya(n+1)))
disp(sprintf('\n   Absolute Relative True Error Percentage'))
disp(sprintf('              = | ( Exact - Approximate ) / Exact | * 100'))
disp(sprintf('              = | %g / %g | * 100',yf-ya(n+1),yf))
disp(sprintf('              = %g',abs( (yf-ya(n+1))/yf )*100))

disp(sprintf('\nThe approximation can be more accurate if we made our step'))
disp(sprintf('size smaller (that is, increasing the number of steps).\n\n'))
