function [t,y]=oderkf(f,t0,y0,tf,ea,er)
%
%        [t,y]=oderkf(f,t0,y0,tf,ea,er)
%
%     Copyright: Eduardo Casas (Octubre, 2002)
%
%  Esta función resuelve el problema de valor inicial
%
%             y'(t)=f(t,y(t))  en [t0,tf]
%             y(t0)=y0
%
%  Para ello se utiliza un método Runge-Kutta-Fehlberg 4(5) o 7(8),
%  según la elección realizada en el programa.
%
%  Los parámetros ea y er son opcionales. Si no se suministran
%  ninguno de los dos valores, se tomará por defecto
%
%             ea=1.d-8 y er=1.d-6
%
%  Si se adjudica ea pero no er, entonces se tomará er=0.
%  Estos valores de ea y er son usados para el control local del
%  error
%             error < tol=ea+er*norm(yn)
%
%  Las variables de salida son:
%
% t   Es el conjunto de instantes tn en los que se ha obtenido las
%     aproximaciones yn de la solución.
%
% y   Es un array que contiene las aproximaciones yn de la solución
%     en los instantes tn.
%
%     La instrucción plot(t,y) permitirá una visualización de la
%     solución.
%
if nargin == 5
   er=0;
elseif nargin == 4
   ea=1.d-8;
   er=1.d-6;
end
rmin=0.5;
rmax=5;
hmin=(tf-t0)/1.d+6;
eps2=eps^2;
%
% Elegir aquí el método de integración junto con la asignación del
% orden p adecuado.
%
[a,b,c,est]=rkf78;
p=0.125;
%[a,b,c,est]=rkf45a;
%[a,b,c,est]=rkf45b;
%p=0.2;
t=t0;
y=y0;
h=max([(ea+er*norm(y0))^p,hmin]);
ind=1;
if t0+h > tf
   h=tf-t0;
   tf=t0+h;
end
while t0 < tf
   [y1,en]=rkf(f,t0,h,y0,a,b,c,est);
   tol=ea+er*norm(y1);
   if en < tol
      ind=ind-1;
      t0=t0+h;
      y0=y1;
      t=[t t0];
      y=[y y0];
      en=max([eps2,en]);
      r=0.9*(tol/en)^p;
      if r < 1
         h=r*h;
      elseif ind <= 0 & r > 1.1
         h=min([rmax r])*h;
      end
      if t0+h > tf
         h=tf-t0;
         tf=t0+h;
      end
   else
      ind=3;
      r=0.9*(tol/en)^p;
      h=min([rmin r])*h;
      if h < hmin
         '       No es posible obtener la precision deseada'
         return
      end
   end
end
%%%%%%%%%%%%%%%%%%%%%%
function [a,b,c,est]=rkf78
c = [0 2/27 1/9 1/6 5/12 1/2 5/6 1/6 2/3 1/3 1 0 1]';
a =[zeros(1,13);2/27 zeros(1,12);1/36 1/12 zeros(1,11);
   1/24 0 1/8 zeros(1,10);5/12 0 -25/16 25/16 zeros(1,9);
   1/20 0 0 1/4 1/5 zeros(1,8);-25/108 0 0 125/108 -65/27 125/54 zeros(1,7);
   31/300 0 0 0 61/225 -2/9 13/900 zeros(1,6);
   2 0 0 -53/6 704/45 -107/9 67/90 3 zeros(1,5);
   -91/108 0 0 23/108 -976/135 311/54 -19/60 17/6 -1/12 zeros(1,4);
   2383/4100 0 0 -341/164 4496/1025 -301/82 2133/4100 45/82 45/164 18/41 0 0 0;
   3/205 0 0 0 0 -6/41 -3/205 -3/41 3/41 6/41 0 0 0;
   -1777/4100 0 0 -341/164 4496/1025 -289/82 2193/4100 51/82 33/164 12/41 0 1 0];
s=41/840;
b =[zeros(1,5) 34/105 9/35 9/35 9/280 9/280 0 s s]';
est=[s zeros(1,9) s -s -s]';
%%%%%%%%%%%%%%%%%%%%%%%%%%
function [a,b,c,est]=rkf45a
c = [0 2/9  1/3  3/4  1  5/6]';
a  = [zeros(1,6);2/9 zeros(1,5);1/12 1/4 zeros(1,4);69/128 -243/128 135/64 0 0 0;
     -17/12 27/4 -27/5 16/15 0 0;65/432 -5/16 13/16 4/27 5/144 0];
b = [47/450 0 12/25 32/225 1/30 6/25]';
est = [1/150 0 -3/100 16/75 1/20 -6/25]';
%%%%%%%%%%%%%%%%%%%%%%%%%%
function [a,b,c,est]=rkf45b
c = [0 1/4  3/8  12/13  1  1/2]';
a = [zeros(1,6);1/4 zeros(1,5);3/32 9/32 zeros(1,4);
    1932/2197 -7200/2197 7296/2197 0 0 0;439/216 -8 3680/513 -845/4104 0 0;
    -8/27 2 -3544/2565 1859/4104 -11/40 0];
b = [16/135 0 6656/12825 28561/56430 -9/50 2/55]';
est = [-1/360 0 128/4275 2197/75240 -1/50 -2/55]';
%%%%%%%%%%%%%%%%%%%%%%%%%%
function [yn1,en]=rkf(f,tn,hn,yn,a,b,c,est)
%
%     Copyright: Eduardo Casas (Octubre, 2000)
%
m=length(b);
kn=feval(f,tn,yn);
for j=2:m
   kn=[kn,feval(f,tn+hn*c(j),yn+hn*(kn*a(j,1:j-1)'))];
end
yn1=yn+hn*(kn*b);
en=hn*norm(kn*est);

