function [tiempo,solE] = Met1Richardson
%
% RESOLUCION PROBLEMA DE CAUCHY SDO (III): METODOS DE 1 PASO +
% EXTRAPOLACION DE RICHARDSON (q = 1/nR)
%
% DATOS A PROPORCIONAR:
% t0:   INSTANTE INICIAL
% T:    LONGITUD DEL INTERVALO
% M:    DIMENSION DEL SISTEMA
% y0:   DATO INICIAL (VECTOR COLUMNA DE DIMENSION M)
% N:    NO. DE PASOS (h = T/N)
% nR:   NO. DE SUBDIVISIONES PARA RICHARDSON (METODO APLICADO CON PASOS h
%       Y qh, CON q = 1/nR; NECESARIAMENTE nR >= 2)
% met:  METODO (1 = EULER + RICHARDSON,
%               2 = RUNGE-KUTTA DE CUARTO ORDEN + RICHARDSON)
% ifig: SALIDA GRAFICA (>  0 GRAFICAS SOLUCIONES CALCULADAS, 
%                            SOLUCION EXACTA, ERROR
%                       >= 0 NO HAY)
%
% ARGUMENTOS DE SALIDA:
% tiempo: INSTANTES DE TIEMPO DONDE SE CALCULA LA SOLUCION 
%                               (TABLERO DE DIMENSION N+1)
% solE:   SOLUCION (TABLERO DE DIMENSION M X (N+1))
%
%
format short;
close all;
%
% DATOS A PROPORCIONAR
%
t0 = 0.;
T = 1.;
M = 1;
%
%y0 = 1.; N = 50;       % TEST 1
%y0 = 1.; N = 50;       % TEST 2
%y0 = 0.; N = 10;       % TEST 3.a
%y0 = 1.e-2; N = 50;    % TEST 3.b
%y0 = 0.; N = 10;       % TEST 4.a
y0 = 0.; N = 50;       % TEST 4.b
%y0 = 0.; N = 100;      % TEST 4.c
%y0 = .2; N = 50;       % TEST 5
%y0 = [.2; 1,]; N = 50; % TEST 6
%
nR = 2;
met = 2;
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(' METODOS POSIBLES: 1 = EULER, 2 = HEUN, 3 = EULER MEJORADO');
%disp('                   4 = RUNGE-KUTTA DE CUARTO ORDEN');
%disp('                   OTRO = RESOLUCION SIMULTANEA');
%met = input(' METODO?  ');
%disp(' ');
%disp(' FIGURAS POSIBLES: 1 = GRAFICA t/y');
%disp('                   2 = CAMPO DE FASES');
%ifig = input(' FIGURAS?  ');
%
% PREPARACION AL CALCULO
%
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);
%
NR = N*nR; q = 1./nR; hR = q*h; nRm1 = nR-1;
hRs2 = .5*hR; hRs6 = hR/6.;
NR1 = NR+1;
qE = q; qRK4 = q^4;
AqE = 1./(1.-q); AqRK4 = 1./(1.-qRK4);
tiempoR = t0:hR:t0+T;
%
sol = zeros(M,N1);
solR = zeros(M,NR1);
solE = zeros(M,N1);
sol(:,1)=y0; solR(:,1)= y0; solE(:,1)= y0;
%
% RESOLUCION
%
if (met==1)
%
% EULER EXPLICITO
%
    iR = 1;
%
    for i=1:N
        i1 = i+1;
        sol(:,i1) = Euler(tiempo(i),h,sol(:,i));
        jR = iR + nRm1;
        for j=iR:jR
            j1 = j+1;
            solR(:,j1) = Euler(tiempoR(j),hR,solR(:,j));
        end;
        solE(:,i1) = AqE*(solR(:,j1) - qE*sol(:,i1));
        iR = iR+nR;
    end;
%
    if (ifig>0)
        solplot = zeros(2,N1);
        for k=1:M
            figure;
            solplot(1,:) = sol(k,:);
            solplot(2,:) = solE(k,:);
            plot(tiempo,solplot,'LineWidth',2);
            h1 = legend('EULER (E)','EULER + RICHARDSON (E+R)',2);
            set(h1,'Interpreter','none');
            xlabel('t'); ylabel('y_k calculada');
        end;
%
% 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
%    solex = ????;                                        % TEST 5    
%    solex = [???? ; ????]                                % TEST 6
        figure;
        plot(tiempo,solex,'LineWidth',2);
        h2 = legend('SOLUCION EXACTA',2);
        set(h2,'Interpreter','none');
        xlabel('t'); ylabel('y EXACTA');
        
        error = zeros(2,N1);
        for k=1:M
            figure;
            error(1,:) = sol(k,:) - solex(k,:);
            error(2,:) = solE(k,:) - solex(k,:);
            plot(tiempo,error,'LineWidth',2);
            h1 = legend('ERROR E','ERROR E+R',2);
            set(h1,'Interpreter','none');
            xlabel('t'); ylabel('error_k');
        end;
    end;
%
elseif (met==2)
%
% RUNGE-KUTTA 4 EXPLICITO
%
    iR = 1;
%
    for i=1:N
        i1 = i+1;
        sol(:,i1) = Runge_Kutta_4(tiempo(i),tiempo(i1),...
            h,hs2,hs6,sol(:,i));
        jR = iR + nRm1;
        for j=iR:jR
            j1 = j+1;
            solR(:,j1) = Runge_Kutta_4(tiempoR(j),tiempoR(j1),...
                hR,hRs2,hRs6,solR(:,j));
        end;
        solE(:,i1) = AqRK4*(solR(:,j1) - qRK4*sol(:,i1));
        iR = iR+nR;
    end;
%
    if (ifig>0)
        solplot = zeros(2,N1);
        for k=1:M
            figure;
            solplot(1,:) = sol(k,:);
            solplot(2,:) = solE(k,:);
            plot(tiempo,solplot,'LineWidth',2);
            h1 = legend('RUNGE-KUTTA 4 (RK4)',...
                'RUNGE-KUTTA 4 + RICHARDSON (RK4+R)',2);
            set(h1,'Interpreter','none');
            xlabel('t'); ylabel('y_k calculada');
        end;
%
% 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
%    solex = ????;                                        % TEST 5    
%    solex = [???? ; ????]                                % TEST 6
        figure;
        plot(tiempo,solex,'LineWidth',2);
        h2 = legend('SOLUCION EXACTA',2);
        set(h2,'Interpreter','none');
        xlabel('t'); ylabel('y EXACTA');
        error = zeros(2,N1);
        for k=1:M
            figure;
            error(1,:) = sol(k,:) - solex(k,:);
            error(2,:) = solE(k,:) - solex(k,:);
            plot(tiempo,error,'LineWidth',2);
            h1 = legend('ERROR RK4','ERROR RK4+R',2);
            set(h1,'Interpreter','none');
            xlabel('t'); ylabel('error_k');
        end;
    end;
%
end;

function yb = Euler(a,h,ya)
yb=ya + h*f(a,ya);

function yb = Euler_mejorado(a,h,hs2,ya)
fa = f(a,ya);
yb = ya + h*f(a + hs2,ya + hs2*fa);

function yb = Heun(a,b,h,hs2,ya)
fa = f(a,ya);
yab = ya + h*fa;
yb = ya + hs2*(fa + f(b,yab));

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



