%% First order differential equations: direction field, numerical solution

clearvars; 

%% You must first download dirfield.m and put this file in the same folder as your mlx-files.
% We consider a first order differential equation of the form
%         y'(t)=f(t,y(t))
% Example:
%        y'(t)=t^2-y(t)^2
% We will usually write this as y'=t^2-y^2.
%% Defining the function f(t,y):
% Here the function f(t,y) is given by
%       f(t,y)=t^2-y^2
% In Matlab we define this function by

f = @(t,y) t^2-y^2;

% Note: we must use f=@(t,y)... even if t or y don't occur in the function.
%% Drawing the direction field
% Consider a solution y(t): At t_0 it has the value y_0=y(t_0), and it must have the slope f(t_0,y_0) at this point.
% We evaluate f(t_0,y_0) on a grid of points in the t,y plane and indicate this slope by a short line:
% For the grid we use the t-values -2,-1.8,...,1.8,2 and the y-values 0, 0.2,...,2.8,3:

dirfield(f,-2:.2:3,0:.2:3)
xlabel('t'); ylabel('y'); title('direction field for f(t,y)=t-y^2')

% Note: Try different values for the step size (here 0.2) to obtain a nice looking direction field.
% A solution of the ODE corresponds to a function where the graph "fits in the direction field".
%
%% Initial value problem
% In order to get a unique solution we must specify an initial condition
%          y(t_0)=y_0
% Here we use the initial condition
%        y(-2)=1
% with t_0=-2 and y_0=1

t0=0; y0=1;
hold on;          % plot over the previous graph
plot(t0,y0,'bo')  % mark initial point (t0,y0) with 'o'

%% Finding an (approximate) solution with ode45/ode15s
%% Solution to the right
% We want to find the solution with the inital t-value t_0 and the final t-value T=3
% Note: we can also have T<t_0, this means we are looking for the solution to the left of the starting point.

T = 3;

% We use the Matlab function ode45:

[ts,ys] = ode45(f,[t0,T],y0);

% This gives a table of values for the solution function y(t):

[ts,ys]          % In Matlab you can click on the array to see all values.
disp([ts,ys])    % this shows the whole array of values (scroll to view all rows)

% The first column contains t-values starting with t_0 and ending with T, the second column contains the corresponding y-values. Use disp to show the whole array (you can scroll to view all values).
% We can plot the solution using

plot(ts,ys,'b')

%
%% Solution to the left
% We now try to find a solution for t between t_0=0 and -2:

[ts,ys] = ode45(f,[t0,-2],y0);

% The Matlab warning tells us that it was not able to find the solution all the way to -2, it was only able to go up to

t1 = ts(end)  % last entry of vector ts
plot(ts,ys,'b')
plot(t1*[1 1],ylim,'k--') % plot vertical line at t=t1 (dashed black)
hold off
title('Solution of IVP y''=t-y^2, y(0)=1')

% Here we show the asymptote at t=t_1 as a dashed vertical line.
% The solution of this IVP exists on the interval (t_1,\infty) with t_1\approx -1.03746
% NOTE: If the solution ends at a point with vertical tangent is is better to use ode15s in place of ode45. In this situation the function ode45 may take a long time and may not recognize that the solution ceases to exist at the point with vertical tangent.
%
%
%
%
%