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='-30*y' ;
   f=inline(fcnstr,'t','y') ;

% t0, initial time

   t0=0 ;

% y0, corresponding value of y at t0

   y0=1.0/3.0;

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

   tf=1.5;

% n, number of steps to take

   n= 10;

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

% 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 ;
Yexact(1) = y0;

for i=1:n
% 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 ;
  Yexact(i+1) = y0*exp(-30*ta(i+1));    % compute exact solution
end


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

% 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(ta,Yexact,'--','LineWidth',2,'Color',[0 0 1]);            
plot(ta,ya,'-','LineWidth',2,'Color',[0 1 0]);
legend('Exact','Approximation');
