%% prog004 metodo de euler

clearvars;
close all;
clc;

%% incializacion

%ecuacion diferencial
% y'' + t*y' + y = 0

%uso de variable auxiliar
% y' = u          <---------------> u = y'
% y'' = -t*y' - y <---------------> u' = t*u  - y


%tiempo incial/final
tini = -0;
tfin =  5;

%valores inciales%
y0 = 1;         %posición incial
u0 = 2;         %velocidad incial u(0) = y'(0)

%incrementos
h1 = 0.5000;
h2 = 0.1000;
h3 = 0.0500;
h4 = 0.0010;

%valores de xn por prelocalización
t1 = tini:h1:tfin;       %valores de t1
t2 = tini:h2:tfin;       %valores de t2
t3 = tini:h3:tfin;       %valores de t3
t4 = tini:h4:tfin;       %valores de t4


N1 = length(t1);     %tamanho de t1
N2 = length(t2);     %tamanho de t2
N3 = length(t3);     %tamanho de t3
N4 = length(t4);     %tamanho de t3

%valores yn por prelocalización
y1 =zeros(1,N1);     %incialmente lleno de ceros
y2 =zeros(1,N2);     %incialmente lleno de ceros
y3 =zeros(1,N3);     %incialmente lleno de ceros
y4 =zeros(1,N4);     %incialmente lleno de ceros

%valores un por prelocalización
u1 =zeros(1,N1);     %incialmente lleno de ceros
u2 =zeros(1,N2);     %incialmente lleno de ceros
u3 =zeros(1,N3);     %incialmente lleno de ceros
u4 =zeros(1,N4);     %incialmente lleno de ceros

%condicion incial de desplazamiento
y1(1) = y0;         %condicion incial desplazamiento
y2(1) = y0;         %condicion incial desplazamiento
y3(1) = y0;         %condicion incial desplazamiento
y4(1) = y0;         %condicion incial desplazamiento

%condicion inicial de velocidad
u1(1) = u0;         %condicion incial velocidad u(0) = y'(0)
u2(1) = u0;         %condicion incial velocidad u(0) = y'(0)
u3(1) = u0;         %condicion incial velocidad u(0) = y'(0)
u4(1) = u0;         %condicion incial velocidad u(0) = y'(0)

%solucion 1
for n = 1:N1-1
    y1(n+1) = y1(n) + h1*u1(n);
    u1(n+1) = u1(n) + h1*( -t1(n)*u1(n) - y1(n) );
end;


%solucion 2
for n = 1:N2-1
    y2(n+1) = y2(n) + h2*u2(n);
    u2(n+1) = u2(n) + h2*( -t2(n)*u2(n) - y2(n) );
end;

%solucion 3
for n = 1:N3-1
    y3(n+1) = y3(n) + h3*u3(n);
    u3(n+1) = u3(n) + h3*( -t3(n)*u3(n) - y3(n) );
end;

%solucion 4
for n = 1:N4-1
    y4(n+1) = y4(n) + h4*u4(n);
    u4(n+1) = u4(n) + h4*( -t4(n)*u4(n) - y4(n) );
end;

figure(1)
plot(t1,y1,'b',t2,y2,'r',t3,y3,'g',t4,y4,'k')
title('solucion euler');
xlabel('t (s)')
ylabel('x(t) (m)');
legend(['y1 = ' num2str(h1)],['y2 = ' num2str(h2)],['y3 = ' num2str(h3)],['y4 = ' num2str(h4)])
axis([tini tfin min([min(y1),min(y2),min(y3),min(y4)]) max([max(y1),max(y2),max(y3),max(y4)]) ]); 
grid on;
box on;
hold on;

%este metodo posee muchas limitaciones inclusive para pequeños valores de h
%puede traer errores hasta del 35% en pocas iteraciones