%% prog004 metodo de euler

%originalmente esta es la ecuación

% dy
%---- + 2*y*t = 0
% dt

% con condicion incial
%y(0) = 1 = exp(0);


%la cual puede ser acomodada

% dy
%---- = f(t,y)= -2*y*t
% dt

clearvars;
close all;
clc;

%% incializacion

%valores inciales
tini = -0;
tfin = 3;
y0 = exp(-tini^2);
t0 = tini;

%incremento
deltat = 1/10;

%valores de tn por prelocalización
tn = tini:deltat:tfin;      %valores de t
Nn = length(tn);            %tamanho de t

%valores yn por prelocalización

yn = zeros(1,Nn);	%incialmente lleno de ceros
yn(1) = y0;         %condicion incial

%% solucion analitica
adeltat = 1/5000;
t_analitic = tini:adeltat:tfin;
y_analitic = exp(-t_analitic.^2);

%% campo direccional y solucion analitica

%espacio de estados inciales
Lt = tfin;
Ly = max(y_analitic);
dt = 0.05;
dy = 0.05;

[t0,y0] = meshgrid(tini:dt:Lt,min(y_analitic):dy:1.1*max(y_analitic));

%velocidad incial
v0 = -2*y0.*t0; %v0 = dy/dt

%vector
figure(1)
quiver(t0,y0,y0,v0,'r');          %dibuja flechas 2D 
title('campo direccional');
xlabel('t_0 (s)');
ylabel('y_0 (m/s)');
legend('campdir')
axis([tini Lt min(y_analitic) 1.1*max(y_analitic)]);
grid on;

figure(2)
quiver(t0,y0,y0,v0,'r');
hold on;
plot(t_analitic,y_analitic,'b');
title('campo direccional y solucion analitica');
xlabel('t_0 (m)');
ylabel('y_0 (m/s)');
legend('campdir','y(t)_{analitic}');
axis([tini Lt min(y_analitic) 1.1*max(y_analitic)]);
grid on;

%% solucion numerica

%solucion 
for n = 1:Nn-1
    fty = -2*yn(n)*tn(n);
    yn(n+1) = yn(n) + fty*deltat;
end;

%% solucion

figure(7)
plot(tn,yn,'k');
hold on;
plot(t_analitic,y_analitic,'b');
title('solucion numerica y solucion analitica');
xlabel('t (s)');
ylabel('y (m)');
legend('y(t)_{numeric}','y(t)_{analitic}');
axis([tini tfin min(yn) 1.1*max(y_analitic)]);
grid on;
box on;

figure(8)
quiver(t0,y0,y0,v0,'r');
hold on;
plot(t_analitic,y_analitic,'b');
hold on;
plot(tn,yn,'k -o');
title('campo direccional solucion analitica y solucion numerica');
xlabel('t_0 (s)');
ylabel('y_0 (m)');
legend('campdir','y(t)_{analitic}','y(t)_{numeric}');
axis([tini Lt min(y_analitic) 1.1*max(y_analitic)]);
grid on;


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