function [tout,yout] = rk23(odefun,tspan,y0,options) % % function [tout,yout] = rk23(odefun,tspan,y0,varargin) % a(2,1) = 1/2; a(3,1) = -1; a(3,2) = 2; c(2) = 1/2; c(3) = 1; b(2) = 1; bhat(1) = 1/6; bhat(2) = 2/3; bhat(3) = 1/6; yout(:,1) = y0; InitialStep = 0.01; RelTol = 1e-3; AbsTol = 1e-6; if (nargin == 4) if (isfield(options,'RelTol')) RelTol = options.RelTol; end if (isfield(options,'AbsTol')) AbsTol = options.AbsTol; end if (isfield(options,'InitialStep')) InitialStep = options.InitialStep; end end k = InitialStep; tout(1) = tspan(1); n = 0; yout(:,n+1) = y0(:); f = zeros(length(y0),3); while (tout(n+1) < tspan(2)) n = n+1; k = min(k,tspan(2)-tout(n)); f(:,1) = odefun(tout(n),yout(:,n)); f(:,2) = odefun(tout(n)+c(2)*k,yout(:,n)+a(2,1)*k*f(:,1)); % rk2 yout(:,n+1) = yout(:,n)+k*(b(1)*f(:,1)+b(2)*f(:,2)); f(:,3) = odefun(tout(n)+c(3)*k,yout(:,n)+a(3,1)*k*f(:,1)+a(3,2)*k*f(:,2)); % prima si calcola la differenza err = k*((bhat(1)-b(1))*f(:,1)+(bhat(2)-b(2))*f(:,2)+bhat(3)*f(:,3)); normyout = norm(yout(:,n+1),Inf); normerr = norm(err,Inf); if (normerr > AbsTol+normyout*RelTol) disp('time step rejected') n = n-1; else % rk3 yout(:,n+1) = yout(:,n+1)+err; tout(n+1) = tout(n)+k; end k = min(2,max(0.6,abs(0.7*... ((AbsTol+normyout*RelTol)/normerr)^(1/3))))*k; end yout = yout.';