function Holmes
% M.H. Holmes, Introduction to Perturbation Methods,
% Springer, New York, 1995, solves the ODE
%
%   y'' + epsilon*lambda*y' + y + ...
%         epsilon*kappa*y^3 = epsilon*cos((1+epsilon*omega)t)
% 
% with parameter values that results in a highly oscillatory
% solution with slowly increasing magnitude.  Here the solution
% is compared to the RMS average computed with ODEAVG.

epsilon = 0.05;
lambda = 2;
kappa = 4;
omega = 3/8;

tspan = [0,250];
y0 = [0;0];

delta = 35;
average = input('Average 1 time or 2 times? ');
avg_opt = [average,2];
sol = odeavg(avg_opt,delta,@ode,tspan,y0);

solode = ode23(@ode,tspan,y0);

% Compute the RMS average.
R = sqrt(sol.y(1,:));
plot(sol.x,R,'r',solode.x,solode.y(1,:),'k')
legend('RMS averaged solution','Computed solution',...
       'Location','SouthWest')
title(['Using \Delta = ',num2str(delta),' and average = ',...
        int2str(sol.average),'.'])

%===Nested function================================================
function yp = ode(t,y)
    yp = [y(2); 0];
    yp(2) = - epsilon*lambda*y(2) - y(1) - epsilon*kappa*y(1)^3 ...
            + epsilon*cos((1+epsilon*omega)*t);
end % ode
%==================================================================

end % Holmes