function test_ode45_1(relerr,abserr)
% Comparación de métodos explícitos para igual evaluaciones de f
 clear all, close all 
 ecdif=@ec04;   % ecuación diferencial a integrar: en ec04.m
 tspan=[0,2];   y0=0;   h=.25;
 % Diversos métodos con misamas llamadaa a f
 [T1,Y1]=eulere(ecdif,tspan,y0,h/4);     % Euler
 [T2,Y2]=heun(ecdif,tspan,y0,h/2);       % Heun
 [T3,Y3]=rk4(ecdif,tspan,y0,h);          % Runge-Kutta orden 4
 options = odeset('RelTol',1.e-9,'AbsTol',1.e-12);
 [T4,Y4]=ode45(ecdif,tspan,y0,options);  % Exacto: R-K de Matlab
 axis([tspan(1), tspan(2), 0, .6]);        % intervalos de ejes
 xlabel('tiempo'); ylabel('y'); plot(T1,Y1,':r');
 hold on; grid
 plot(T2,Y2,'--g');   plot(T3,Y3,'-.b');   plot(T4,Y4,'-k'); 
 title('Solución de y''=y(1/t-2t)'); grid on
 legend('Euler','Heun','rk4','ode45');
end

function yprima=ec04(t,y)
% solución exacta desconocida
 if t==0, yprima=1;  else  yprima=y*(-2*t+1/t);  end
end

function [T,Y]=eulere(f,tspan,y0,h)  % Método de Euler
 n=fix((tspan(2)-tspan(1))/h+1);
 T=tspan(1)*zeros(n,1); Y=y0*zeros(n,length(y0));
 y=y0; t=tspan(1);
 for i=2:n
    y1=y+h*f(t,y);  t1=t+h;   % Fórmula de Euler
    T(i,:)=t1; Y(i,:)=y1';    % tiempo valor de la funcion
    t=t1; y=y1;               % se actualizan t e y
 end
end

function [T,Y]=heun(f,tspan,y0,h)    % Método de Heun
 n=(tspan(2)-tspan(1))/h+1;
 T=tspan(1)*zeros(n,1); Y=y0*zeros(n,length(y0));
 y=y0; t=tspan(1);
 for i=2:n
    t1=t+h;
    yp=f(t,y);                % predictor con Euler explícito
    y1=y+h*yp;
    y1n=y+h*(yp+f(t1,y1))/2;  % corrector con la regla trapezoidal
    T(i,:)=t1; Y(i,:)=y1n;    % tiempo valor de la funcion
    t=t1; y=y1;               % se actualizan t e y
 end
end

function [T,Y]=rk4(f,tspan,y0,h)     % Runge-Kutta orden 4
 n=fix((tspan(2)-tspan(1))/h+1);
 T=tspan(1)*zeros(n,1); Y=y0*zeros(n,length(y0));
 y=y0; t=tspan(1);
 for i=2:n
    t1=t+h;
    k1=h*f(t,y);              % cálculo de los coeficientes k
    k2=h*f(t+h/2,y+k1/2);
    k3=h*f(t+h/2,y+k2/2);
    k4=h*f(t+h,y+k3);
    y1=y+(k1+2*k2+2*k3+k4)/6;
    T(i,:)=t1; Y(i,:)=y1;     % tiempo valor de la funcion
    t=t1; y=y1;               % se actualizan t e y
 end
end