function [t,yx,y,y_h,y_hi] = Heun_4(dydt,tspan,y0,h,es,maxit)
% Método de Heun para integrar una ODE
%  dydt: función que evalua la ODE; tspan: intervalos de cálculo
%  y0: punto de partida; h: paso; es: tolerancia de iteración corrector
%  maxit: máximo de iteraciones; yx: valor exacto; y: valor con Euler
%  y_h: valor con Heun sin itera.; y_hi: con iteración
 if nargin<6, maxit = 50; end; if nargin<5, es = 0.0001; end 
 ti = tspan(1); tf = tspan(2); t = (ti:h:tf)'; n = length(t);
 if t(n)<tf, t(n+1)=tf; n=n+1; end                  % Por si acaso, un intervalo más en t
 y = y0*ones(n,1); y_h = y; y_hi = y;               % Inicializaar valores
 yx = 4/1.3*(exp(0.8*t)-exp(-0.5*t))+2*exp(-0.5*t); % Valor exacto
 for i = 1:n-1
    hh = t(i+1) - t(i); iter = 0; 
    k1 = dydt(t(i),y(i)); k2 = dydt(t(i),y_h(i)); k3 = dydt(t(i),y_hi(i)); % Derivadas
    y(i+1)    = y(i)   + k1*hh;                                            % Euler
    y_h(i+1)  = y_h(i) + hh/2*(k2+dydt(t(i+1),y_h(i)+hh*k2));     % Heun sin iteración
    y_hi(i+1) = y_hi(i)+ k3*hh;
    while 1
       yold= y_hi(i+1); k2 = dydt(t(i+1),y_hi(i+1));
       y_hi(i+1) = y_hi(i) + (k3+k2)*hh/2;                        % Heun con iteración
       iter = iter + 1;
       if y_hi(i+1)~=0, ea=abs((y_hi(i+1)-yold)/y_hi(i+1))*100; end
       if ea <= es || iter >= maxit, break, end
    end                       
 end
 plot(t,yx,'b--',t,y,'r-',t,y_h,'k--',t,y_hi,'g:'); xlabel('t'); ylabel('y');
 title('Solución de y''=4e^{(0.8t)}-0.5y'); grid on
 legend('Exacto','Euler','Heun','Heun con it.','Location','NorthWest');
end

function dp = dpdt_Heun_c(t, p)
 dp = 4*exp(0.8*t)-0.5*p;
end