%Euler
clear all
%Datos
%Introducimos los valores de las constantes
a=0.003;
b=0.3;
c=0.01;
h=input('Introducir valores de h:'); 
t0=0;
y0=input('Introducir vector [So,Io]:');
tN=40;
% calculamos los subintervalos
t=t0:h:tN;          % con linspae si nos diesen Nº de elementos
N=(tN-t0)/h;
y=zeros(2,length(t));
y(:,1)=y0';
for i= 1:N;
    y(:,i+1)=y(:,i)+h*[-a*y(1,i).*y(2,i);a*y(1,i).*y(2,i)-b*y(2,i)-c*y(2,i)];
end

% grafico
plot(t,y)
I=y(2,:);
[fila,col]=find(I==max(max(I)))
Diademaximo=(col-1)*h
legend('Euler','Location','best'); %Lo ultimo es para que la legenda no interfiera
