%% Eulerexample 2
%%
% 
% $$y'= -ty+1/y.^2, y(1)=1$$
% on $[1,2]$
% 
clear; close all;format compact
%% 
f=@(t,y) -t.*y+1./y.^2
help myEuler
%%
Trange=[1 2];y0=1;
n=10;
[T,Y]=myEuler(f,Trange,y0,n);
plot(T,Y,'-*');grid on
title(['Euler''s method, y''-ty+1/y.^2, y(1)=1)']); shg
hold on
%% Vary intial value
% Repeat
n=2*n
figure 
hold on
for y0=0.1:.1:1
    [T,Y]=myEuler(f,Trange,y0,n);
    plot(T,Y,'-',T,Y,'.')
end    
%% More ... repeat this line
figure
n=2*n;[T,Y]=myEuler(f,Trange,0.1,n); plot(T,Y,'.--');
figure
[T,Y]=myEuler(f,[1 1.1],0.1,n); plot(T,Y,'.--');
h=.1/n
%% Compare to ode45
Trange=[1 2];y0=0.1;
[T,Y]=myEuler(f,Trange,y0,n);
figure
plot(T,Y,'-',T,Y,'.');grid on;shg
%% Change Trange:
Trange=[0 1];y0=1;
n=10;
[T,Y]=myEuler(f,Trange,y0,n);
plot(T,Y,'-*');grid on
title(['Euler''s method, y''-ty+1/y.^2, y(0)=1)']); shg
hold on

%% Symbolic: Leave this!
syms y(t)
Dequ=diff(y)==t*y+1/y.^2
analy=dsolve(Dequ(t),y(1)==1);
%ysym=vectorize(analy)
%%
%t=linspace(1,2);
%Y=((exp((3*t.^2)/2).*(2*exp(-3/2) + 6^(1/2)*pi^(1/2)*erf((6^(1/2).*t)/2)...
% - 6^(1/2)*pi^(1/2)*erf(6^(1/2)/2)))/2)^(1/3)
% syms vähän tökkii
%plot(t,Y)
%
%% Change into
%
% $$y'= -ty-1/y.^2, y(1)=1$$
%
%%
clear close all
f=@(t,y) -t.*y-1./y.^2
Trange=[1 2];y0=1;
n=10;
[T,Y]=myEuler(f,Trange,y0,n);
plot(T,Y,'-*');grid on
title(['Euler''s method, y''-ty-1/y.^2, y(1)=1)']); shg
hold on
%%
n=2*n,[T,Y]=myEuler(f,Trange,y0,n);plot(T,Y,'-*');shg
% Blows up at about 1.23
%% What does ode45 say?
figure
[T,Y]=ode45(f,Trange,y0);
plot(T,Y,'.',T,Y);grid on;shg
%{
Warning: Failure at t=1.226114e+00.  Unable to meet integration tolerances without reducing the step size below
the smallest value allowed (3.552714e-15) at time t.
%}
%% 
%% Direction field
close all
n=16;
tpoints=linspace(1,2,n); 
% ypoints=linspace(.1,2,n);  % Too much variation in lengths.
ypoints=linspace(.5,2,n);
[t,y]=meshgrid(tpoints,ypoints); % Rectangular grid, remember 3d-graphics.
%%
plot(t(:),y(:),'.');shg
%% 
hold on
f=@(t,y) -t.*y-1./y.^2 
vt=ones(size(y));
vy=f(t,y);   % The derivative of solution curve, remember: y'=f(t,y)
hold on
quiver(t,y,vt,vy,1.5);shg
%% This doesn't help in the different sizes case
%
norms=sqrt(1+vy.^2);
quiver(t,y,vt/norms,vy/norms,5);shg

%%

title('y''= t+y','FontSize',14)
%title('y''=-y-5e.^{-t}sin(5t)','FontSize',14)
xlabel('t');ylabel('y','Rotation',0)
xlim([0 3.2]);ylim([-1.3 1.15]); % Tune axis limits


Yfinal=[ya(end);yb(end);yc(end)]
Tfinal=ta(end)   % ==tb(end) == tc(end)
yexactFinal=exp(Tfinal)-Tfinal-1 
errors=yexactFinal-[ya(end);yb(end);yc(end)]
diff(errors)
