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

This is an autonomous ODE

with the function

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

for

. Hence we get stationary points

(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

. Let us plot the solutions for

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

solving for y0=-0.5:

Warning: Failure at t=-1.098608e+00.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t.

solving for y0=0:
solving for y0=0.5:
solving for y0=1:
solving for y0=1.5:

Warning: Failure at t=1.098610e+00.  Unable to meet integration tolerances without reducing the step size below the smallest value allowed (3.552714e-15) at time t.

hold off; title('solutions for five initial conditions y(0) = -.5, 0, .5, 1, 1.5')

We see that

the solution for blows up at , it exists for
the solution for blows up at , it exists for
the other three solutions exist on the whole real axis

We have the initial condition

. For a given T we want to find

. Integrating

from

to T gives

Example: Let

. We will get

. We can find this T by computing the improper integral

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
T =           1.09861228866811

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