function mnhf_rk4_euler(h,Nmax)
%MNHF_RK4_EULER illustrates use of RK4 in solving the non-autonomous ODE 
% y'=1+sin(\pi*y)+exp(-t) with y(0)=1. 
% Solution obtained by RK4 is contrasted with that obtained using explicit 
% Euler.
%
% h - time step
% Nmax - number of iterations
%
% Usage  mnhf_rk4_euler(0.1,100)

%% Initialize independent variable array.
t = linspace(0.0,h*Nmax,Nmax+1);

%% Initialize dependent variable arrays; specify initial condition.
yRK4 = zeros(size(t)); yRK4(1) = 1.0;
yEE = zeros(size(t));  yEE(1) = 1.0;

for ij=1:Nmax
 %% Calculate solution based on explicit Euler.
 yEE(ij+1) = yEE(ij)+h*(1.0+sin(pi*yEE(ij))+exp(-t(ij)));
 %% Calculate solution based on RK4.
 k1 = 1.0+sin(pi*yRK4(ij))+exp(-t(ij));
 k2 = 1+sin(pi*(yRK4(ij)+0.5*h*k1))+exp(-(t(ij)+0.5*h));
 k3 = 1+sin(pi*(yRK4(ij)+0.5*h*k2))+exp(-(t(ij)+0.5*h));
 k4 = 1+sin(pi*(yRK4(ij)+h*k3))+exp(-(t(ij)+h));
 yRK4(ij+1) = yRK4(ij)+h/6.0*(k1+2.0*k2+2.0*k3+k4);
end

%% Plot results.
figure(1); hold on; box on
set(gca,'fontsize',18)
plot(t,yEE,'k--','linewidth',2)
plot(t,yRK4,'k-','linewidth',2)
xlabel('{\it t}'); ylabel('{\it y}')
legend('explicit Euler','RK4','location','best')
title(['time step, {\it h} = ' num2str(h)])