%% Simulation 2 Companion
%% MATLAB Onramp
% In order to understand Simulation 2, you should first complete the following modules from the MATLAB Onramp.
% Programming
% It will be assumed that you are familiar with the content of these modules in the discussion below.
%
%% Slope Fields
% Let us use MATLAB to plot slope fields.
% Goal: Plot the slope field for \frac{dy}{dt} = t(3-y).
% First I provide you with a function that produces a slope field. You are free to try to understand the makeup of this function using the provided comments but your task will only require you call the function.
% Note: You should not uncomment anything in this box as what is needed has already been uncommented.

%==========================================================================
%
% slopefield: Draws the slope field of a first-order ordinary differential equation dy/dt = f(t,y).
%
%    slopefield(f,tval,yval)
%
%--------------------------------------------------------------------------
%
% ------
% INPUT:
% ------
%   f : (function_handle) dy/dt = f(t,y)
%   tval : (row vector) of EVENLY-SPACED t-values at which to sketch slopes
%   yval: (row vector) of EVENLY–SPACEd y-values at which to sketch slopes
%
%==========================================================================

function slopefield(f,tval,yval)

% calculates increments in t and y assuming the increments in t are equal
% and the increments in y are equal
tinc = tval(2)-tval(1);
yinc = yval(2)-yval(1);

% length of slope lines to ensure no overlap 
slopelength = min(tinc,yinc)/2;

% creates grid (T,Y) of coordinates
[T,Y]=meshgrid(tval,yval); 

% evaluates the slopes dy/dt = f(T,Y)
dydt = f(T,Y); 

% identifies components of slope lines 
dT = slopelength*cos(atan(dydt));
dY = slopelength*sin(atan(dydt));

% draws thick red slope lines from (T-dT/2,Y-dT/2) to (T+dT/2,y+dT/2)
quiver(T,Y,dT/2,dY/2,0,'.r',LineWidth=1); 
hold on
quiver(T,Y,-dT/2,-dY/2,0,'.r',LineWidth=1);
hold off

% sets axes limits
Tmin = tval(1)-1/2*tinc;
Tmax = tval(end)+1/2*tinc;
Ymin = yval(1)-1/2*yinc;
Ymax = yval(end)+1/2*yinc;
axis([Tmin Tmax Ymin Ymax])

end

% Next we create a function handle. You write:
% function_name = @(input_variable1,input_variable2,...) (function_formula)
% Reminder: uncomment each line start and select Run Section.

% f = @(t,y) (t.*(3-y)) % define f(t,y) = t(3-y)

% A function handle is convenient, for example because it can be used in normal function notation later. For example look at the output for the following.

% f(1,2) % should yield f(1,2) = 1(3-2) = 1

% Next I want to define the time interval to be from 0 to 10 with 1 spacing and the y interval to be from 0 to 6 with 1 spacing.

% tval = 0:1:10
% yval = 0:1:6

% Now we call the slopefield function that had been defined earlier using our function handle f and our equally spaced intervals t and y.

% slopefield(f,tval,yval) % slopefield for dy/dt = t(3-y)

% Now ain't that a beauty.
%% Solution Curves
% Now suppose we want to sketch a solution curve given an initial value.
% Goal: Plot solution curves to \frac{dy}{dt} = t(3-y) at various initial values.
% We will work with ode45 which is a numerical solver based on an the Runge-Kutte method. You should not uncomment anything in this box.

%==========================================================================
%
% ode45: Numerically solves the initial value problem dy/dt = f(t,y) with
% y(0) = y0.
%
%    [t,y] = ode45(f,tspan,y0)
%
%--------------------------------------------------------------------------
%
% ------
% INPUT:
% ------
%   f : (function_handle) dy/dt = f(t,y)
%   tspan = [t0 tf] numerically solves from t0 to tf
%   y0 : (number/row-vector) initial value/s
%
% ------
% OUTPUT:
% ------
%   t : (column vector) t-values in numerical solution 
%   y : (column vector) values y(t) in numerical solution
%
%==========================================================================

% Let us continue with the function handle f for f(t,y) = t(3-y) we have defined earlier and find the solution to dy/dt = f(t,y) with initial value y(0) = 1 with t in [0,10].

% y0 = 1 % initial value
% tspan = [0, 10] % domain
% [t,y] = ode45(f,tspan,y0) % solve: you can see the output vectors t and y on the right

% Next let us plot.

% hold on % plot on the slopefield we had created earlier
% plot(t,y,Linewidth=1) % plot with Linewidth=1
% hold off

% What if I want to plot various the solution for various initial values? Let use use for to plot the solutions for all initial 0,1,...,6. To understand the for loop it was expected you completed the Programming Matlab Onramp as indicated at the top of the page.

% for y0 = 0:6
%     [t,y] = ode45(f,tspan,y0); % solve the `;` hides the output
%     hold on % plot on the slopefield we had created earlier
%     plot(t,y,Linewidth=1) % plot with Linewidth=1
% end
% hold off

% Wow that is beautiful. Look at how all the solutions converge to the solution y = 3 which is what we would call an equilibrium solution.