%% Runge-Kutta Method
% Solves the initial value problem y'= x + y, y(x0)=y0 

clc, clear all, close all
%% 
% Define the RHS of the ODE
f=@(x,y) x+y; 
% Initial values
x0=0; y0=1;
% Final value of x
xf=1;
% Step size, 
h=0.05; 
% number of steps
N=floor((xf-x0)/h);

%%
x=x0; y=y0;
X=x0; Y=y0; % stores the initial values
hold on

for i=1:N
    k1 = f(x,y);
    k2 = f(x+h/2,y+h/2*k1);
    k3 = f(x+h/2,y+h/2*k2);
    k4 = f(x+h,y+h*k3);
    
    y = y + h/6*(k1+2*k2+2*k3+k4);
    x = x + h;
  
    plot(x,y,'o') % plots the computed approximation
    
    X=[X;x]; Y=[Y;y]; % stores all computed values 
end
grid on
xlabel('x'); ylabel('y')
title('Runge-Kutta for the solution of dy/dx=x+y')

%% 
% If an explicit solution can be found, it can eb entered below
% for comparing it with the approximate solution, otherwise
% the entire code below must be commented out! 
yexact = @(x) -x-1+2*exp(x);
% Plot of the exact solution for comparison purposes only
fplot(yexact, [x0 xf],'Color','r')
hold off

