function [t,y]=eulerrichardson(f,a,b,ya,nsteps)
% Euler's method for ODE IVP's, dy/dt=f(t,y) y(a)=ya, a <= t <= b,
% with Richardson extapolation.
% Author: Brenton LeMesurier
% Date: March 8, 2007
% Version for MATLAB7 only: function input as a function pointer.
%
% Input Variables
% - f is (a pointer to) the function f(t,y).
% - a, b specify the initial and final t values.
% - nsteps is the number of time steps to take,
% so step size is dt=(b-a)/nsteps.
%
% Output Variables
% - t is an array ot time values, a:dt:b
% - y is an array of corresponding values of the solution y(t)

% For a Matlab 6 version: input f as a formula in a text string instead, and use this inline commmand.
% f=inline(f,'t','y');
% Create an arrays t and y. These could be created on the fly, but doing it
% this way shows the form of the arrays, and is slightly faster.
[thalfsteps,ymoreaccurate]=euler(f,a,b,ya,2*nsteps);
[t,ylessaccurate]=euler(f,a,b,ya,nsteps);
% extract the elements of the better (smaller dt) solution that are for
% times also present in the less accurate one:
ymoreaccuratesome=ymoreaccurate(1:2:end); % "end" means the last passobe value of the index
y=2*ymoreaccuratesome-ylessaccurate; % Euel is frst order accurate, so "2-1", not "4-1" as for second order.