Code in Matlab: Complete lines s1,s2, s3, p(p=3) , yn1 and err so that the ode324 code implements the embedded explicit trapezium and simpson pair.

Code:

% INPUTS: A function handle f = @(t,y) ... such that dy/dt = f(t,y) % The times span tspan(1) < t

% Initial step size h = sqrt(tol); % Initialise solver t = tspan(1); tend = tspan(2); yn = y0; flag = 0; % Storage tt = t; yy = yn;

%% Loop through time while ( t < tend ) % RK stages s1 = ; s2 = ; s3 = ; p = ; yn1 = yn + ; % Error estimate err = ; % Step size adjust: hnew = 0.8*(tol*min(abs(yn1)./err))^(1/(p+1))*h; hnew = min(hnew, 2*h); % Don't let h grow too fast. if ( norm(err, inf) > tol ) % Error is too large: if ( ~flag ) % Try the new time step, flag = 1; % Set flag to remember failure. else % Time step failed. hnew = h/2; % Take half of old one, end else % Error is small enough! % Store solution and move to next time: yn = yn1; t = t + h; yy = [yy, yn1]; tt = [tt, t]; flag = 0; % Successful step. Set flag to happy. end % Choose next step (ensure we finish at the end of the interval!) h = min(hnew, tend - t);

end

if ( nargout == 0 ) subplot(2,1,1) plot(tt, yy(1,:), '.'), shg title('First component of solution') subplot(2,1,2) semilogy(tt(2:end), abs(diff(tt))) title(['No. of time steps = ' int2str(length(tt))]) tt = []; yy = []; end

end