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.
%Here f = (f1, f2)^T, y = (y1, y2)^T
% The code is for solving Illustration problem of Section 5.11

%  define f1(t,y)
   fcnstr='9*y1 + 24*y2 + 5*cos(t) - 1/3*sin(t)' ;
   f1=inline(fcnstr,'t','y1', 'y2') ;
   
%  define f2(t,y)
   fcnstr2='-24*y1 -51*y2 -9*cos(t) + 1/3*sin(t)' ;
   f2=inline(fcnstr2,'t','y1', 'y2') ;
   
% t0, initial time
   t0=0;

% y1_0, corresponding value of y1 at t0
   y1_0=4/3;

% y2_0, corresponding value of y1 at t0
   y2_0=2/3;

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

   tf=1.0;

% n, number of steps to take

   n=10;

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

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


ta(1)=t0 ;
ya1(1)=y1_0 ;
ya2(1)=y2_0 ;

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

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

% Calculating k1_1, k2_1, k3_1, and k4_1 for y1, and 
% Calculating k1_2, k2_2, k3_2, and k4_2 for y1

  k1_1 = f1(ta(i),ya1(i), ya2(i)) ;
  k1_2 = f2(ta(i),ya1(i), ya2(i)) ;
  
  k2_1 = f1(ta(i)+0.5*h,ya1(i)+0.5*k1_1*h,  ya2(i)+0.5*k1_2*h) ;
  k2_2 = f2(ta(i)+0.5*h,ya1(i)+0.5*k1_1*h,  ya2(i)+0.5*k1_2*h) ;
  
  k3_1 = f1(ta(i)+0.5*h,ya1(i)+0.5*k2_1*h, ya2(i)+0.5*k2_2*h) ;
  k3_2 = f2(ta(i)+0.5*h,ya1(i)+0.5*k2_1*h, ya2(i)+0.5*k2_2*h) ;
  
  k4_1 = f1(ta(i)+h,ya1(i)+k3_1*h,  ya2(i)+k3_2*h) ;
  k4_2 = f2(ta(i)+h,ya1(i)+k3_1*h,  ya2(i)+k3_2*h) ;
     
% Using 4th Order Runge-Kutta formula
  ya1(i+1)=ya1(i)+1/6*(k1_1+2*k2_1+2*k3_1+k4_1)*h ;
  ya2(i+1)=ya2(i)+1/6*(k1_2+2*k2_2+2*k3_2+k4_2)*h ;

   
end

% compute exact solution at these time points
for i=1:n+1
   y1(i)=2.0*exp(-3*ta(i)) - exp(-39*ta(i)) + 1/3*cos(ta(i));
   y2(i)=-exp(-3*ta(i)) + 2.0* exp(-39*ta(i)) - 1/3*cos(ta(i));
end 
  
% Plotting the Exact and Approximate solution of the ODE.
hold on
xlabel('t');ylabel('y');
title('Exact and Approximate Solution of the ODE by the 4th Order Runge-Kutta Method');
plot(ta,y1,'--','LineWidth',2,'Color',[0 0 1]);            
plot(ta,ya1,'-','LineWidth',2,'Color',[0 1 0]);
plot(ta,y2,'--','LineWidth',1,'Color',[0 1 1]);            
plot(ta,ya2,'-','LineWidth',1,'Color',[0 1 0.5]);

legend('Exact-1','Approximation-1', 'Exact-2','Approximation-2');