%% tarea01 EDO05 metodo de euler

%originalmente esta es la ecuación

% dy
%---- - ty^2 + y/t = 0
% dt

% con condicion incial
%y(1) = 1


%la cual puede ser acomodada

% dy
%---- = f(t,y)= ty^2 - y/t
% dt

clearvars;
close all;
clc;

%% incializacion

%valores inciales
tini = 1;
tfin = 2;
y0 = 1;
t0 = tini;

%incremento
deltat = 1/100;

%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 numerica

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


%% campo direccional y solucion analitica

%espacio de estados inciales
Lt = tfin;
Ly = max(yn);
dt = 0.01;
dy = 0.01;

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

%velocidad incial
v0 = t0.*(y0.^2) - y0./t0;%  ty^2 - y/t; %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(yn)/1.1 1.1*max(yn)]);
grid on;



%% solucion

figure(2)
plot(tn,yn,'-o b','Linewidth',1);
title('solucion numerica');
xlabel('t (s)');
ylabel('y (m)');
legend('y(t)_{numeric}');
axis([tini tfin min(yn)/1.1 1.1*max(yn)]);
grid on;
box on;

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


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