function [T,Y]=EM(u0,a,f,k)
% u0 = condicion inicial;
% a = alpha
% f = numero real de la EDO
% k = numero k, con el que se obtiene el paso h( h=2^k)
% L = hasta donde se aproxima la EDO
L=800;
% c = constante de la EDO
c=0.175;
% h = paso de la aproximacion
h=2^k;
% T = vector que contiene los t con que se evalua la
%funcion solucion de la Edo
T=0:h:L;
%Y= vector que contiene los valores aproximados
% de la EDO evaluados en x
Y=zeros(1,length(T));
% Se añade condicion inicial al vector Y
Y(1)=u0;
% Se establece u0 como primer valor de y
y=u0;
% Se itera hasta llegar a L
for n=1:length(T)-1
    % Se calcula el valor del integrando en el punto medio
    % del intervalo [n*h,(n+1)*h]
    Yy=y+(h/2)*(c*(y-a)*(1-y)+f);
    % Se calcula y(n+1) valor de la aproximacion de la EDO
    % usando Yy
    y=y+h*(c*(Yy-a)*(1-Yy)+f);
    % Se agrega aproximacion n+1-esimo de y al vector Y
    Y(n+1)=y;
end
end

------------------------------------------------------------------------------------------------------------------------------

function e=Error(k,Y)
% Se Hace una aproximacion de la EDO con el metodo
% de Euler Modificado y se guarda la aproximacion
% en el vector y
[t,y]=EM(0.2,0.1,0,k);
% Se almacena el maximo error en la variable e
e=0;
% Se busca el maximo |yref(tn)-y(n)|, para tn=hn<=800
for i=1:length(t)-1
    % Se calcula el yref(tn)
    Yref=Y(i*(2^(10+k)));
    % Se guarda el posible nuevo maximo en A
    A=abs(Yref-y(i));
    % Se verifica que A es el nuevo maximo
    if A>e
        % Se reemplaza e por el valor de A
        e=A;
    end
end


----------------------------------------------------------------------------------------------------------------

% Se ejecuta el algoritmo Euler Modificado para
% un paso h=2^(-10) y se guarda en T e Y la
% aproximacion obtenida
[T,Y]=EM(0.2,0.1,0,-10);
% Crear un vector de errores para los nueve h(pasos)
E=zeros(1,9);
% Crear un vector para los nueve h (pasos)
H=zeros(1,9);
% Se calcula el error para los nueve pasos
for k=-5:3
    % Se guardan los pasos en H
    H(k+6)=2^k;
    % Se guardan los errores en E
    E(k+6)=Error(k,Y);
end
%empezar a graficar
hold on
%graficar
plot(H,E,'redo')
% Título del grafico
title('Gráfico de error e(h) versus paso h');
% Label del eje x
xlabel('Paso del método de Euler Modificado [u]');
% Label del eje y
ylabel('Error del método de Euler Modificado');
% Terminar de graficar
hold off
% Estimacion del orden del metodo en funcion de
% e(h, p guarda la suma para calcular el promedio
% del orden
p=0;
for k=2:9
    %suma a p la i-esima estimacion del orden
    p=p+(log(E(k)/E(k-1)))/(log(2));
end
%calcula el promedio de los ordenes
p=1.0*p/8.0


-----------------------------------------------------------------------------------------------

function f=minf(C,k)
% c = numero al que se quiere que converja la funcion
% k = numero que determina el paso h (h=2^k)
% e = error con el que se evaluara la convergencia
e=0.0001;
% ad = paso con el que se busca f
ad=0.0001;
% Lo siguiente hace una linea en y=1 para
% mostrar que la aproximacion no converge
w=1:800;
z=ones(1,length(w));
plot(w,z,'green-')
% h= paso del metodo de Euler Modificado
h=(2^k);
% Se elige f inicial 0.0175 pues desde ese punto la EDO
% aproximada no se va a -infinito
f=0;
% Elegimos un un valor de f donde se detendra la busqueda
% para que no sea un loop infinito la busqueda en caso que
% no exista f tal que la Edo converge. Elegimos 4 pues si
% no encontramos f antes tampoco lo encontraremos despues
% pues con 4 la EDO converge a más de 3
while f<=4
    hold on
    % Se hace una aproximacion a la EDO para verificar si con f
    % la aproximacion converge a C, guardamos la aproximacion en
    % Y
    [T,Y]=EM(0,0.1,f,k);
    plot(T,Y,'-')
    xlim([0 800]);
    ylim([-2 2]);
    % La revision de "convergencia" parte de n=1
    n=1;
    % Se utilizan n's hasta el que indexa al penultimo elemento de
    % Y
    while h*n<800
        % Se verifica que el n-esimo termino de Y este "cerca" de
        % C y del siguiente termino
        if (C-e)<=Y(n) && Y(n)<=(C+e) && (Y(n)-e)<=Y(n+1)&& Y(n+1)<=(Y(n)+e)
            % Se verifica que cada elemento que sigue del que cumplió
            % la condicion anterior este "cerca" de C
            while (C-e)<=Y(n) && Y(n)<=(C+e)
                % si hasta el ultimo termino de Y se culple lo anterior
                % se concluye que la aproximacion converge a C
                if n*h==800
                    return
                end
                % Se aumenta el indice n para iterar
                n=n+1;
            end
        end
        % Se aumenta el indice n para iterar
        n=n+1;
    end
    % Si la aproximacion no "converge" buscamos un f aumentado
    % en ad y se verifica si converge a C
    f=f+ad
end
end

------------------------------------------------------------------------------------------------------------
function H=maxh
% h = el paso con el cual se hara la aproximacion
% a la EDO. Se iniciará con h=0.00001
h=0.0001;
% e = error con el cual se verificara que la
% aproximacion "converge" a 1 y no diverge. Se determinó
% que tiende a 1 pues al ejecutar EM con diferentes
% pasos la solucion se acerca a 1
e=3.0001;
% C = valor al que se acercará la aproximacion
C=0;
% ad = paso con el que se busca h
ad=0.001;
% a = constante de la EDO
a=0.1;
% u0 = condicion inicial de la EDO
u0=2*a;
% f = constante de la EDO
f=0;
% Se itera infinitamente hasta encontrar el h.
% Se podria poner un limite a la iteracion, pero
% como así arroja un h es innecesario
while true
    % Dado que a la funcion EM se le ingresa una
    % variable para obtener el paso (h=2^k) se
    % debe despejar k
    k=log(h)/log(2);
    % Se hace una aproximacion a la EDO para verificar si con h
    % la aproximacion es estable, guardamos la aproximacion en
    % Y
    [T,Y]=EM(u0,a,f,k);
    %  La revision de "estabilidad" parte de n=1
    n=1;
    % Se itera la revision de "estabilidad" con los indices de de Y
    while n<=length(T)-1
        % Si el ultimo elemento de Y no esta cerca de C entonces
        % H es el maximo paso donde la aproximacion es estable
        if (Y(n)<(C-e) || (C+e)<Y(n)) && n==length(T)-1
            return
        end
        % Se verifica que cada elemento Y(n), con n entre 1 y length(T)-1
        % este "cerca" de C, si es así entonces la aproximacion esta es
        % estable
        while (C-e)<=Y(n) && Y(n)<=(C+e) && n<=length(T)-1
            % Si se pasa la prueba anterior, hasta n=length(T)-1, la
            % aproximacion es "estable" con h y se tiene que h es el
            % nuevo maximo paso tal que la aproximacion es estable
            if n==length(T)-1
                % Se guarda h como nuevo posible maximo-h-estable
                H=h
                break
            end
            % Se aumenta el indice n en 1 para iterar
            n=n+1;
        end
        % Se aumenta el indie n en 1 para iterar
        n=n+1;
    end
    % Si h es estable, h se aumenta en ad para buscar donde
    % h deja de ser estable
    h=ad+h;
end