function test_ode113(lambda,relerr,abserr)

% This is a test program for the ode solver 'ode113'.
% The test is carried out for the single equation
%     y' = lambda*y + (1-lambda)*cos(t) - (1+lambda)*sin(t)
% The initial value at t=0 is y(0)=1.  The true solution is
%     y = cos(t) + sin(t)
% The user can input the relative and absolute error
% tolerances to be used by ode113.  These are incorporated
% using the initialization program 'odeset'. 
% The program can be adapted easily to other equations and 
% other parameter values.

% Initialize and solve
options = odeset('RelTol',relerr,'AbsTol',abserr);
t_begin = 0; t_end = 20;
y_initial = true_soln(t_begin);
num_fcn_eval = 0;  % initialize count of derivative evaluations
soln = ode113(@deriv,[t_begin,t_end],y_initial,options);...
                              % See below for function deriv.

% Produce the solution on a uniform grid using interpolation
% of the solution obtained by ode113. The points plotted with
% 'o' are for the node points returned by ode113.
h_plot = (t_end-t_begin)/1000; t_plot = t_begin:h_plot:t_end;
y_plot = deval(soln,t_plot);
figure
plot(soln.x,soln.y,'o',t_plot,y_plot)
title(['Interpolated solution:',...
    ' points noted by ''o'' are at ode113 solution nodes'])
xlabel(['\lambda  = ',num2str(lambda)])
disp('press on any key to continue')
pause 

% Produce the error in the solution on the uniform grid.
% The points plotted with 'o' are for the solution values
% at the points returned by ode113.
y_true = true_soln(t_plot);
error = y_true - y_plot;
y_true_nodes = true_soln(soln.x);
error_nodes = y_true_nodes - soln.y;
figure
plot(soln.x,error_nodes,'o',t_plot,error)
title('Error in interpolated solution')
xlabel(['\lambda  = ',num2str(lambda)])

norm_error = norm(error,inf);
disp(['maximum of error = ',num2str(norm_error)])
disp(['number of derivative evaluations = ',...
    num2str(num_fcn_eval)])

function dy = deriv(t,y)
% ----------------------

% Define the derivative in the differential equation.
dy = lambda*y + (1-lambda)*cos(t) - (1+lambda)*sin(t);
num_fcn_eval = num_fcn_eval + 1;

end % deriv

function true = true_soln(t)
% ------------------------

% Define the true solution of the initial value problem.
true = sin(t) + cos(t);

end % true_soln

end % test_ode113