% This file executes the Runge-Kutta method for a 1st-order ODE


clear all

h = .2;                 % Choose the step size

x0 = 0;                 % Choose initial x value
xf = 1;                % Choose the final x value. (Note that xf - x0 must
                        % be divisible by h in order to end exactly at xf)

X = x0:h:xf;            % This sets equally spaced x values
N = length(X);          % We will use the length of the vector of x values
Y = zeros(size(X));     % Preallocates space for y values


Y(1) = 1;               % Choose initial y value


for j = 1:N-1
    k1 = frhs(X(j),Y(j));                   % left hand slope f(X(j),Y(j))
    k2 = frhs(X(j)+.5*h,Y(j)+.5*h*k1);      % first middle slope f(X(j)+h/2,Y(j) + h/2*k1)
    k3 = frhs(X(j)+.5*h,Y(j)+.5*h*k2);      % second middle slope f(X(j)+h/2,Y(j) + h/2*k2)
    k4 = frhs(X(j)+h,Y(j)+h*k3);            % right hand slope f(X(j+1),Y(j) + h*k3) 
    Y(j+1) = Y(j) + h/6*(k1 + 2*k2 + 2*k3 + k4);  % This is the next step approximation
end



Yexact = exp(X);    %because our DE is y'(x)=y, with y(0)=1

figure(1)
plot(X,Yexact, X,Y)
legend('Exact solution', 'Approximated solution')



error = abs(Yexact(N) - Y(N))
relerror = error/abs(Yexact(N))