function [tiempo,sol] = Met1SDO
%
% RESOLUCION PROBLEMA DE CAUCHY SDO (I): METODOS DE 1 PASO
%
% 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)
% met: METODO (1 = EULER,
%              2 = EULER MEJORADO,
%              3 = HEUN,
%              4 = RUNGE-KUTTA DE CUARTO ORDEN)
%              OTRO = RESOLUCION SIMULTANEA Y ANALISIS DEL ERROR)
% ifig: SALIDA GRAFICA (1 = GRAFICA SOLUCION,
%                       2 = CAMPO DE FASES CUANDO M = 2)
%
% 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.;
M = 1;
%
y0 = 1.; N = 50;       % TEST 1
%y0 = 1.; N = 50;       % TEST 2
%y0 = 0.; N = 50;       % 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
%
met = 5;
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);
format short;
close all;
%
% RESOLUCION
%
if (met==1)                                                                % EULER EXPLICITO
%
% EULER EXPLICITO
%
    sol = zeros(M,N1);
    sol(:,1)=y0;
%
    for i=1:N
        sol(:,i+1) = Euler(tiempo(i),h,sol(:,i));
    end;
%
    figure;
    if (ifig==1)
        plot(tiempo,sol);                                                  % FIGURA SOLUCION
    elseif (ifig==2)
        plot(sol(1,:),sol(2,:))                                            % FIGURA CAMPO DE FASES
    end;
%
    legend('EULER EXPLICITO',2);
%
elseif (met==2)                                                            % EULER MEJORADO
%
% EULER MEJORADO
%
    sol = zeros(M,N1);
    sol(:,1)=y0;
%
    for i=1:N
        sol(:,i+1) = Euler_mejorado(tiempo(i),h,hs2,sol(:,i));
    end;
%
    figure;
    if (ifig==1)
        plot(tiempo,sol);                                                  % FIGURA SOLUCION
    elseif (ifig==2)
        plot(sol(1,:),sol(2,:))                                            % FIGURA CAMPO DE FASES
    end;
%
    legend('EULER MEJORADO',2);
%
elseif (met==3)                                                            % HEUN (EULER MODIFICADO)
%
% HEUN
%
    sol = zeros(M,N1);
    sol(:,1)=y0;
%
    for i=1:N
        sol(:,i+1)=Heun(tiempo(i),tiempo(i+1),h,hs2,sol(:,i));
    end;
%
%
    figure;
    if (ifig==1) 
        plot(tiempo,sol);                                                  % FIGURA SOLUCION
    elseif (ifig==2)
        plot(sol(1,:),sol(2,:))                                            % FIGURA CAMPO DE FASES
    end;
%
    legend('HEUN',2);
%
elseif (met==4)                                                            % RUNGE-KUTTA 4
%
% RUNGE-KUTTA 4
%
    sol = zeros(M,N1);
    sol(:,1)=y0;
%
    for i=1:N
        sol(:,i+1) = Runge_Kutta_4(tiempo(i),tiempo(i+1),...
            h,hs2,hs6,sol(:,i));
    end;
%
    figure;
    if (ifig==1)
        plot(tiempo,sol);                                                  % FIGURA SOLUCION
    elseif (ifig==2)
        plot(sol(1,:),sol(2,:))                                            % FIGURA CAMPO DE FASES
    end;
%
    legend('RUNGE KUTTA 4',2);
%
else
%
% RESOLUCION SIMULTANEA
%
    sol1 = zeros(M,N1);
    sol2 = zeros(M,N1);
    sol3 = zeros(M,N1);
    sol4 = zeros(M,N1);
    sol1(M,1)=y0; sol2(M,1)=y0; sol3(M,1)=y0; sol4(M,1)=y0;
%
    for i=1:N
        sol1(:,i+1) = Euler(tiempo(i),h,sol1(:,i));
        sol2(:,i+1) = Euler_mejorado(tiempo(i),h,hs2,sol2(:,i));
        sol3(:,i+1) = Heun(tiempo(i),tiempo(i+1),h,hs2,sol3(:,i));
        sol4(:,i+1) = Runge_Kutta_4(tiempo(i),tiempo(i+1),...
            h,hs2,hs6,sol4(:,i));
    end;
%
    figure;
%
    solplot = zeros(4,N1);
    solplot(1,:) = sol1(1,:);
    solplot(2,:) = sol2(1,:);
    solplot(3,:) = sol3(1,:);
    solplot(4,:) = sol4(1,:);
    plot(tiempo,solplot);
    h1 = legend('EULER','EULER MEJORADO','HEUN','RUNGE-KUTTA 4',2);
    set(h1,'Interpreter','none');
    xlabel('t'); ylabel('y_1 calculada');
%
    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
%    solex = ????;                                        % 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 = zeros(4,N1);
    for i=1:4
        error(i,:) = solplot(i,:) - solex;
    end;
%
    plot(tiempo,error);
    h3 = legend('ERROR EULER','ERROR EULER MEJORADO','ERROR HEUN',...
        'ERROR RUNGE-KUTTA 4',3);  
    set(h3,'Interpreter','none');
    xlabel('t'); ylabel('error y_1');
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



