% This file executes the Improved Euler's 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             % for loop, with index j iterating from 1 to N-1
    k1 = frhs(X(j),Y(j));        % The RHS of the DE, f(X(j),Y(j))
    k2 = frhs(X(j)+h,Y(j)+h*k1);     % right hand slope estimate
    k = (k1+k2)/2;         % improved euler estimate for slope
    Y(j+1) = Y(j) + h*k;    % next y-value.
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))