%% Direction field and ODE-solver
% 
close all
n=16;
tpoints=linspace(0,3,n); ypoints=linspace(-1,1,n);
[t,y]=meshgrid(tpoints,ypoints); % Pistehila, muistele 3d-grafiikkaa.
%%
plot(t(:),y(:),'.')
%%
pt=ones(size(y));
py=-y-5*exp(-t).*sin(5*t);
hold on
quiver(t,y,pt,py,1.5);
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
%% Solve the equation y'=f(t,y)=y-5*exp(-t)sin(5*t)
f=@(t,y) -y-5*exp(-t).*sin(5*t)
[T,Y]=ode45(f,[0,3],0);
plot(T,Y)

%%
% help ginput
ty = [];  % Gather chosen initial points
n = 0;
% Loop, picking up the points.
disp('Move cursor into graphics window')
disp('******************************')
disp(' ')
disp('Left mouse button picks points.')
disp('Right mouse button picks last point.')
but = 1;
while but == 1
    [ti,yi,but] = ginput(1);
    [T,Y]=ode45(f,[ti,3],yi);
    plot(T,Y)
    n = n+1;
    ty(:,n) = [ti;yi];
end
plot(ty(1,:),ty(2,:),'b*')  % Add chosen initial points to fig.
