%% Direction field and ODE-solver, ODEdirfield.m
% 
% 
%%
% $$y'= f(t,y), f(t,y) = t+y$$
% 
% $$y'= f(t,y), f(t,y) = -y-5e^{-t}sin(5t)$$
% 
%%
close all; clear
%%
n=16;
tpoints=linspace(0,3,n); ypoints=linspace(-1,1,n);
[t,y]=meshgrid(tpoints,ypoints); % Rectangular grid, remember 3d-graphics.
%%
plot(t(:),y(:),'.');shg
%% Define your differential equation here. 
%  
f=@(t,y) t+y         % y' = t+y
% f=@(t,y) -y-5*exp(-t).*sin(5*t)    

%%
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);
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
%% Solve the equation y'=f(t,y)=t+y
f=@(t,y) t+y  % Defined already, just repeat to be sure.
[T,Y]=ode45(f,[0,3],0);
plot(T,Y,'k')  
% Repeat this line on command window changing initial value (possibly color):
[T,Y]=ode45(f,[0,3],1); plot(T,Y,'k')   % Use up-arrow.
% Or CTR-ENTER here
% Use for-loop, perhaps.
%%%% Interactively choose initial values on graphics window. ginput
%  help ginput
t= [];  % 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);
  [Tr,Yr]=ode45(f,[ti,3],yi); % t= ti...3 (up to right window border)
    plot(Tr,Yr)
  [Tl,Yl]=ode45(f,[ti 0],yi); % t = ti...0 (down to left win border)
    plot(Tl,Yl)
    n = n+1;
    ty(:,n) = [ti;yi];
end
plot(ty(1,:),ty(2,:),'b*')  % Add chosen initial points to fig.

%% Solve the equation y'=f(t,y)=y-5*exp(-t)sin(5*t)
f=@(t,y) -y-5*exp(-t).*sin(5*t)  % Defined already above.
[T,Y]=ode45(f,[0,3],0);
plot(T,Y,'k')
%%  More initial values:
[T,Y]=ode45(f,[0,3],1);
plot(T,Y,'k')
%%
% Repeat this line on command window changing initial value (possibly color):
[T,Y]=ode45(f,[0,3],1); plot(T,Y,'k')   % Use up-arrow.
% Or CTR-ENTER here
%% Now, choose initial values pointing on graphics window.
% 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);
[Tr,Yr]=ode45(f,[ti,3],yi); % t= ti...3 (up to right window border)
    plot(Tr,Yr)
  [Tl,Yl]=ode45(f,[ti 0],yi); % t = ti...0 (down to left win border)
    plot(Tl,Yl)
    n = n+1;
    ty(:,n) = [ti;yi];
end
plot(ty(1,:),ty(2,:),'b*')  % Add chosen initial points to fig.
