function [tiempo,sol] = MetAdams
%
% RESOLUCION PROBLEMA DE CAUCHY SDO (II): METODOS DE ADAMS, 2 PASOS (r = 2)
% INICIALIZACIÓN CON RUNGE-KUTTA 4
%
% DATOS A PROPORCIONAR:
% t0: INSTANTE INICIAL
% T: LONGITUD DEL INTERVALO
% y0: DATO INICIAL (VECTOR COLUMNA DE DIMENSION M)
% M: DIMENSION DEL SISTEMA
% N: NO. DE PASOS (h = T/N)
% met: METODO  1 = ADAMS-BASHFORD (AB)
%              2 = ADAMS-MOULTON (AM)
%              3 = PREDICTOR-CORRECTOR ADAMS-BASHFORD-MOULTON (ABM)
% epsm: TOLERANCIA TEST DE PARADA (AM, ABM) 
% iammax: MAXIMO DE ITERACIONES (AM, ABM)
% ifig:  SALIDA GRAFICA  0 = NO
%                        1 = GRAFICA SOLUCION
% idisp: IMPRESION ITERACIONES  0 = NO
%                               1 = ITER + ERROR
%
% ARGUMENTOS DE SALIDA:
% tiempo: INSTANTES DE TIEMPO DONDE SE CALCULA LA SOLUCION 
%                               (TABLERO DE DIMENSION N+1)
% sol: SOLUCION  TABLERO DE DIMENSION M X (N+1) SI met <= 3;
%                TABLERO DE DIMENSION M X 4 X (N+1) SI met = 4
%

%
% DATOS A PROPORCIONAR
%
t0 = 0.;
T = 1.;
%
%y0 = 1.; M = 1; N = 50;         % TEST 1
%y0 = 1.; M = 1; N = 50;         % TEST 2
%y0 = 0.; M = 1; N = 50;         % TEST 3.a
%y0 = 1.e-2; M = 1; N = 50;      % TEST 3.b
%y0 = 0.; M = 1; N = 10;         % TEST 4.a
%y0 = 0.; M = 1; N = 50;         % TEST 4.b
%y0 = 0.; M = 1; N = 100;        % TEST 4.c
y0 = .2; M = 1; N = 50;         % TEST 5
%y0 = [.2; 1,]; M = 2; N = 50;   % TEST 6
%
met = 3;
epsm = 1.e-5; iammax = 100;      % PRECISION Y MAXIMO DE ITERACIONES (AM)
idisp = 1;
ifig = 1;
%
%disp(' ');
%t0 = input(' INSTANTE INICIAL?  ');
%disp(' ');
%T = input(' LONGITUD DEL INTERVALO TEMPORAL?  ');
%disp(' ');
%M = input(' DIMENSION DEL SISTEMA?  ');
%disp(' ');
%y0 = input(' DATO INICIAL? (VECTOR COLUMNA)  ');
%disp(' ');
%N = input(' NO. DE ETAPAS DE TIEMPO?');
%disp(' ');
%disp(' METODO: 1 = ADAMS-BASHFORD,');
%disp('         2 = ADAMS-MOULTON,');
%disp('         3 = PREDICTOR-CORRECTOR ADAMS-BASHFORD-MOULTON');
%met = input(' METODO?  ');
%disp(' ');
%disp(' FIGURAS POSIBLES: 0 = SIN GRAFICA);
%disp('                   1 = GRAFICA t/y');
%disp('                   2 = CAMPO DE FASES');
%ifig = input(' FIGURAS?  ');

%
% PREPARACION AL CALCULO
%
format short;
close all;
%
h = T/N;                                         % PASO DE TIEMPO
hs2 = .5*h; hs6 = h/6.;
N1 = N+1;
tiempo = t0:h:t0+T;                              % PARTICION DEL INTERVALO
%
disp(' TIEMPOS:');
disp(tiempo);
%
bcoef0 = 23./12.; bcoef1 = -4./3.; bcoef2 = 5./12.; % COEFICIENTES AB
bcoef0h = h*bcoef0; bcoef1h = h*bcoef1; bcoef2h = h*bcoef2;
mcoef0 = 19./24.; mcoef1 = -5./24.; mcoef2 = 1./24.;% COEFICIENTES AM
mcoefm1 = 3./8.;                                    % COEFICIENTES AM
mcoef0h = h*mcoef0; mcoef1h = h*mcoef1; mcoef2h = h*mcoef2;
mcoefm1h = h*mcoefm1;

%
% INICIALIZACION: RUNGE-KUTTA 4
%
sol = zeros(M,N1);
sol(:,1) = y0;
sol(:,2) = Runge_Kutta_4(tiempo(1),tiempo(2),h,hs2,hs6,sol(:,1));
sol(:,3) = Runge_Kutta_4(tiempo(2),tiempo(3),h,hs2,hs6,sol(:,2));
fnm2 = f(tiempo(1),sol(:,1));
fnm1 = f(tiempo(2),sol(:,2));
yn = sol(:,3); fnm0 = f(tiempo(3),yn);

%
% RESOLUCION
%
if (met==1)                                    % ADAMS-BASHFORD, r = 2
    for i=4:N1
        sol(:,i) = yn + bcoef0h*fnm0 + bcoef1h*fnm1 + bcoef2h*fnm2;
        fnm2 = fnm1; fnm1 = fnm0;
        yn = sol(:,i); fnm0 = f(tiempo(i),yn);
    end;
%
    if (ifig==1)
        figure
        plot(tiempo,sol);                      % FIGURA SOLUCION
        legend('ADAMS-BASHFORD',2);
    end;
%    
elseif (met==2)                                % ADAMS-MOULTON, r = 2
    for i=4:N1
        errorm = 1.;
        zn1 = yn;
        ti = tiempo(i);
%
        iam = 0;
        while (errorm > epsm)
            fn1 = f(ti,zn1);
            yn1 = yn + mcoefm1h*fn1+ mcoef0h*fnm0 + mcoef1h*fnm1...
                + mcoef2h*fnm2;
            errorm = norm(yn1-zn1); zn1 = yn1;
            iam = iam+1;
            if (iam >= iammax)
                disp(' ADAMS-MOULTON - DEMASIADAS ITERACIONES - ITER =');
                disp(iam);
                break;
            end
        end
%
        if (idisp > 0)
            disp(' ADAMS-MOULTON - ITER, ERROR =');
            disp([iam, errorm]);
        end
%
        sol(:,i) = yn1;
        fnm2 = fnm1; fnm1 = fnm0;
        yn = yn1; fnm0 = f(ti,yn);
    end
%
    if (ifig==1)
        figure
        plot(tiempo,sol);                      % FIGURA SOLUCION
        legend('ADAMS-MOULTON',2);
    end;
%    
elseif (met==3)                                % PREDICTOR-CORRECTOR, r = 2
    for i=4:N1
        errorm = 1.;
        zn1 = yn + bcoef0h*fnm0 + bcoef1h*fnm1 + bcoef2h*fnm2;
        ti = tiempo(i);
%
        iam = 0;
        while (errorm > epsm)
            fn1 = f(ti,zn1);
            yn1 = yn + mcoefm1h*fn1+ mcoef0h*fnm0 + mcoef1h*fnm1...
                + mcoef2h*fnm2;
            errorm = norm(yn1-zn1); zn1 = yn1;
            iam = iam+1;
            if (iam >= iammax)
                disp(' ADAMS-MOULTON, DEMASIADAS ITERACIONES. iam =');
                disp(iam);
                break;
            end
        end
%
        if (idisp > 0)
            disp(' ADAMS-MOULTON - ITER, ERROR =');
            disp([iam, errorm]);
        end
%
        sol(:,i) = yn1;
        fnm2 = fnm1; fnm1 = fnm0;
      
        yn = yn1; fnm0 = f(ti,yn);
    end
%
    if (ifig==1)
        figure
        plot(tiempo,sol);                      % FIGURA SOLUCION
        legend('ADAMS-BASHFORD-MOULTON',2);
    end;
end

%
% EVENTUALMENTE: FIGURAS ADICIONALES
%
if (ifig==1)
    figure;
%
% SOLUCION EXACTA
%
%    c = 1.; solex = exp(c*tiempo);                       % TEST 1
%    solex = -exp(-tiempo)+tiempo.*tiempo-2.*tiempo+2;    % TEST 2
%    solex = .5*tiempo.*tiempo-tiempo;                    % TEST 3.a    
%    solex = 1.e-2*exp(10.*tiempo) + .5*tiempo.*tiempo...
%        - tiempo;                                        % TEST 3.b    
%    aux1 = 1.e4/(1.e4 + 1.); aux2 = 1.e2/(1.e4 + 1.); 
%    solex = aux1*sin(tiempo)-aux2*cos(tiempo)...
%            +aux2*exp(-100.*tiempo);                     % TESTS 4.a,b,c
    r=8.; K=8./13.; Ky0=K*y0; 
    solex = (Ky0*exp(r.*tiempo))./(K+y0*(exp(r.*tiempo)-1.)); % TEST 5    
%    solex = [???? ; ????]                                % TEST 6
%
    plot(tiempo,solex);
    h2 = legend('SOLUCION EXACTA',2);
    set(h2,'Interpreter','none');
    xlabel('t'); ylabel('y_1 exacta');
%
    figure;
%
    error = sol - solex;
    plot(tiempo,error);
    h3 = legend('ERROR',3);
    set(h3,'Interpreter','none');
    xlabel('t'); ylabel('error');
end;

function yb=Runge_Kutta_4(a,b,h,hs2,hs6,ya)
c = a + hs2;
k1=f(a,ya);
k2=f(c,ya+hs2*k1);
k3=f(c,ya+hs2*k2);
k4=f(b,ya+h*k3);
yb=ya + hs6*(k1+2*k2+2*k3+k4);

function z = f(t,y) 
%
%c=1; z=c*exp(c*t); % sol: exp(c*t) para y(0) = 1. % TEST 1
%z=t*t-y;          % sol: -exp(-t)+t.*t-2*t+2  para y(0) = 1. % TEST 2
%z=10.*y+11.*t-5.*t.*t-1.; % sol: t.*t/2-t para y(0) = 0. 
                        % sol: eps*exp(10.*t)+.5*t.*t-t para y(0) = eps
                        % ATENCION: MAL CONDICIONAMIENTO!
                        % TESTS 3.a, 3.b
%z=100.*(sin(t)-y); % sol: 1.e4/(1.e4 + 1)*sin(t)
                   % -1.e2/(1.e4 + 1)*cos(t)
                   % + 1.e2/(1.e4 + 1)*exp(-100.*t)
                   % para y(0) = 0.
                   % ATENCION: VER COMPORTAMIENTO PARA N = 10, 50, 100
                   % TESTS 4.a, 4.b, 4.c
r=8.; K=8./13.; z=r*y*(1-y/K); % sol: ????      % TEST 5
%b=.35; z(1)=y(2); z(2)=-2.*b*y(2)-y(1)+cos(t);  % 
                   % sol: ????                   % TEST 6



