% =====================================================================================================
% =====            Discritized Approximation of The Exact Solution of: d^2y/dx^2 = x+y            =====
% =====       The Exact Solution of The Ordinary Differential Equation is: (e^x)-(e^(-x))-x       =====
% ===== << The Runge - Kutta Method for Higher Than First Order Ordinary Differential Equation >> =====
% =====================================================================================================

clc
clear

dt=0.1;			%The Time Step
Time=0:dt:1;	%The Time Interval

FirstOrderDerivative(1)=1;	%The Initial Condition for y'(0)=1
SecondOrderDerivative(1)=0; %The Initial Condition for y(0)=0
Z(1)=FirstOrderDerivative(1);

for i=1:length(Time)-1
   
   K1=dt*Z(i);
   X(i)=SecondOrderDerivative(i);
   
   TheFunction1(i)=SecondOrderDerivative(i)+X(i);
   L1=dt*TheFunction1(i);
   
   TheFunction2(i)=Z(i)+((1/2)*L1);
   K2=dt*TheFunction2(i);
   
   TheFunction3(i)=(X(i)+((1/2)*dt))+(SecondOrderDerivative(i)+((1/2)*K1));
   L2=dt*TheFunction3(i);
   
   TheFunction4(i)=Z(i)+((1/2)*L2);
   K3=dt*TheFunction4(i);
   
   TheFunction5(i)=(X(i)+((1/2)*dt))+(SecondOrderDerivative(i)+((1/2)*K2));
   L3=dt*TheFunction5(i);
   
   TheFunction6(i)=Z(i)+L3;
   K4=dt*TheFunction6(i);
   
   TheFunction7(i)=(X(i)+dt)+(SecondOrderDerivative(i)+K3);
   L4=dt*TheFunction7(i);
   
   SecondOrderDerivative(i+1)=SecondOrderDerivative(i)+((1/6)*(K1+2*K2+2*K3+K4));
   
   Z(i+1)=Z(i)+((1/6)*(L1+2*L2+2*L3+L4));

end

ExactSolution=exp(Time)-exp(-Time)-Time; %The Exact Solution

plot(Time,ExactSolution,'b',Time,SecondOrderDerivative,'r');
legend('The Exact Solution','The Approximate Solution');
title('Higher Runge - Kutta Approximation Method, dt=0.1');
xlabel('Time');
ylabel('Y');
grid on