function [tiempo,sol] = MetML
%
% RESOLUCION PROBLEMA DE CAUCHY SDO (II): METODO MATLAB
%
%   ODE23 is an implementation of the explicit Runge-Kutta (2,3) pair of
%   Bogacki and Shampine called BS23.  It uses a "free" interpolant of order
%   3.  Local extrapolation is done.

%   Details are to be found in The MATLAB ODE Suite, L. F. Shampine and
%   M. W. Reichelt, SIAM Journal on Scientific Computing, 18-1, 1997.

%
% 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.;
tfinal = t0 + T;
M = 1;
%
%y0 = 1.;       % TEST 1
%y0 = 1.;       % TEST 2
%y0 = 0.;       % TEST 3.a
y0 = 1.e-2;    % TEST 3.b
%y0 = 0.;       % TEST 4.a
%y0 = 0.;       % TEST 4.b
%y0 = 0.;       % TEST 4.c
%y0 = .2;       % TEST 5
%y0 = [.2; 1,]; % TEST 6
%
iplot = 1; idisp = 0;
%
%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
%
format short;
close all;
% options = odeset('RelTol',1.e-3,'Abstol',1.e-6); VALORES POR DEFECTO
options = odeset('RelTol',1.e-6,'Abstol',1.e-9);
%
% RESOLUCION
%
[tiempo,sol] = ode23(@frate,[t0,tfinal],y0,options);
%
if (iplot > 0)
    figure;
    plot(tiempo,sol);                                                      % FIGURA SOLUCION
    h1 = legend('ODE 23 - MATLAB',2);
    set(h1,'Interpreter','none');
    xlabel('t'); ylabel('y 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 exacta');
%
    figure;
%
    error = sol - solex;
    plot(tiempo,error);
    h3 = legend('ERROR',2);  
    set(h3,'Interpreter','none');
    xlabel('t'); ylabel('error');
end
%
if (idisp > 0)
    disp(' '); disp(' TIEMPOS:'); disp(tiempo);
    disp(' '); disp(' SOLUCION:'); disp(sol);
end;

function [z] = frate(t,y)
%z = zeros(size(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



