%% Example for 1st order autonomous ODE
% Here we use numerical Matlab commands.

clearvars; clf                                % clear variables, clear graphics
warning('off','MATLAB:fplot:NotVectorized')   % get rid of stupid warning
format long g                                 % show 15 digits of numerical results

% You first have to download the file dirfield.m and put it in the same folder as your mlx files.
% We consider the ODE
%         y'(t) = y(t)\cdot \bigl(y(t)-1\bigr)
% This is an autonomous ODE y'(t)=g(y(t)) with the function g(y)=y(1-y)

g = @(y) y*(y-1);
fplot(g,[-.5,1.5]); 
set(gca,'XAxisLoc','origin','YAxisLoc','origin'); box off  % draw axes thru origin
xlabel('y'); title('function g(y) = y\cdot(1-y)')

% We see that g(y)=0 for y=0,\; y=1. Hence we get stationary points y=0,\;y=1 (constant solutions).
% We plot the direction field:

f = @(t,y) g(y);
dirfield(f,-2:.25:2,-2:.2:3); hold on
xlabel('t'); ylabel('y'); title('direction field for y'' = y\cdot(1-y)')

%
% We have the initial condition y(0)=y_0. Let us plot the solutions for y_0=-0.5,\;0,\;0.5,\;1,\;1.5:

for y0 = -.5:.5:1.5
  plot(0,y0,'bo')                      % mark initial point with 'o'
  fprintf('solving for y0=%g:\n',y0)
  [ts,ys] = ode45(f,[0,2],y0);         % solve for t going from 0 to 2
  plot(ts,ys,'b')
  if ts(end)<2                         % if solution does not exist up to 2
    plot(ts(end)*[1 1],ylim,'k--')     %  plot vertical asymptote
  end
  [ts,ys] = ode45(f,[0,-2],y0);        % solve for t going from 0 to -2
  plot(ts,ys,'b')
  if ts(end)>-2                        % if solution does not exist up to 2
    plot(ts(end)*[1 1],ylim,'k--')     %  plot vertical asymptote
  end
end
hold off; title('solutions for five initial conditions y(0) = -.5, 0, .5, 1, 1.5')

% We see that
% the solution for y_0=1.5 blows up at a_*\approx 1.0986, it exists for t\in (-\infty,a_*)
% the solution for y_0=-0.5 blows up at b_*\approx -1.0986, it exists for t\in (b_*,\infty)
% the other three solutions exist on the whole real axis
%
% We have the initial condition y(0)=y_0. For a given T we want to find Y=y(T). Integrating \frac{y'(t)}{g(y(t))}=1 from t=0 to T gives
%           \displaystyle
% \int_{y_0}^Y \frac1{g(y)}dy = T
% Example: Let y_0=1.5. We will get y(t)\to\infty as t\to T. We can find this T by computing the improper integral
%           \displaystyle
% \int_{y_0}^\infty \frac1{g(y)}dy = T
% Note that this improper integral is finite, hence the solution blows up at a finite time T.

v = @(y) 1/g(y);
T = integral(v,1.5,Inf,'Arr',true)  % compute integral, use Inf for infinity

% ('Arr',true  allows us to define the function with * / ^  instead of .* ./ .^)