%% MATLAB help file: S19 version
% This file provides examples of some of the MATHLAB commands that we use
% in our math 216 labs.  These are listed by the lab in which they
% first appear.  
%
% You can search this file to find them as well, of course.  in
% alphabetical order, they are @, :, axis, dirfieldsys, disp, eig,
% eulermethod_w19, hold, legend, ode45, ode15s, ones, plot, plot3,
% PlotComet_3D, tic (and toc), title, xlabel, ylabel
%
%-------------------------------------------------------------------------------
%% Lab 0
%----------------------------------------
% plot( xvector, yvector )
%    plot the values in yvector against the corresponding values in
%    xvector
%    example command
xvector = 0:.05:10;
yvector = sin(xvector) ./ xvector;
plot( xvector, yvector );

% to change plot style
plot( xvector, yvector, '-b', 'LineWidth', 2 );

%----------------------------------------
% hold on
% hold off
%    add subsequent plot commands to an existing figure, or stop
%    doing so
%    example command
hold on;
plot( xvector, sin(xvector)./(xvector.^2), '--k', 'LineWidth', 2 );

%----------------------------------------
% title('title text')
% xlabel('label text')
% ylabel('label text')
%    add a title and axis labels to an existing figure
%    example command
xlabel( 'x' );
ylabel( 'y' );
title( 'Plot of sin(x)/x, sin(x)/x^2' );

%----------------------------------------
% legend('label1', 'label2',... )
%    add a legend to an existing figure
%    example command
legend( 'sin(x)/x', 'sin(x)/x^2' );

% if we plot something after generating the legend, it will be 
%    omitted from the legend
plot( [ 0 10 ], [ 0 0 ], ':k' );

%----------------------------------------
% axis( [ xmin xmax ymin ymax ] )
%    reset the axis scales so that xmin <= x <= xmax, ymin <= y <= ymax
%    example command
axis( [ 0 10 -1, 5 ] );

%----------------------------------------
% min:step:max
%    generate a vector [ min min+step min+2*step ... max ]
%    example command
t = 0:0.01:10;

%----------------------------------------
% [ tvals, yvals ] = ode45( function_handle, t_range, initial_condition )
%    generate a numerical approximation to a differential equation
%    example command
[ tvals, yvals ] = ode45( @(t,y) sin(y)/(t+1), [0 10], 1 );

%----------------------------------------
% f = @(variables) expression
%    generate a function handle 
%    example command
f = @(t, y) sin(y)/(t+1);
f(1,1)

%-------------------------------------------------------------------------------
%% Lab 1
%----------------------------------------
% in Lab 1 we use the same commands as we had in Lab 0.

%-------------------------------------------------------------------------------
%% Lab 2
%----------------------------------------
% dirfieldsys( f_handle, g_handle, [xmin,xmax], [ymin,ymax] )
%     example command
dirfieldsys2( @(x,y) y, @(x,y) x.^2 - x.*y, [-3,3], [-3,3] );

%----------------------------------------
% [ tvals, yvals ] = ode45( function_handle, t_range, initial_condition )
%    example command, to solve a system
[ tvals, xvals ] = ode45( @(t,x) [ x(2); x(1)-4*x(2) ], [0 10], [1; 0] );

% ode15s has the same arguments.  In most cases you want ode45, not this
[ tvals, xvals ] = ode15s( @(t,x) [ x(2); x(1)-4*x(2) ], [0 10], [1; 0] );

%----------------------------------------
% plot( xvector, yvector )

% plot a phase portrait for a system given the solution xvals from
% ode45 in a new figure
figure();
plot( xvals(:,1), xvals(:,2), '-k', 'LineWidth', 2 );

%-------------------------------------------------------------------------------
%% Lab 3
%----------------------------------------
% [ evecs, evals ] = eig( matrix )
% evals = eig( matrix )
%     example commands
[ evecs, evals ] = eig( [ 2 -3; 3 2 ] )
evals = eig( [ 1 -3; 3 1 ] )

%----------------------------------------
% format long
%     sets display of output to show more digits
%     example command
format long

%-------------------------------------------------------------------------------
%% Lab 4
%----------------------------------------
% disp( 'text' )
%     example command
disp( 'This text appears in the Command Window.' );

%----------------------------------------
% [tvals, xvals] = eulermethod_w19( function_handle, trange, init, step)
%     example commands
[ tvals, yvals ] = eulermethod_w19( @(t,y) sin(y)/(t+1), [0 10], [1; 0], 0.1 );
[ tvals, xvals ] = eulermethod_w19( @(t,x) [x(2); -x(1)], [0 10], [1; 0], 0.1 );

%----------------------------------------
% [ tvals, yvals ] = ode45( function_handle, t_range, initial_condition )
%    example command
[ tvals, yvals ] = ode45( @(t,y) sin(y)/(t+1), [0 10], 1 );
%    to solve a system
[ tvals, xvals ] = ode45( @(t,x) [ x(2); x(1) 4*x(2) ], [0 10], [1; 0] );

% we can set options as well
% opt = odeset( 'optionName', value,... );
opt1 = odeset( 'MaxStep', 0.1 );
opt2 = odeset( 'RelTol', 1e-8 );
opt3 = odeset( 'MaxStep', 0.1, 'RelTol', 1e-8 );
[tLS, xLS] = ode45(...
    @(t,x) [10*(-x(1)+x(2)); 28*x(1)-x(2)-x(1)*x(3); -(8/3)*x(3)+x(1)*x(2)],...
    [0 2000], [5;5;5], opt2 );

% ode15s has the same arguments.  In most cases you want ode45, not this
[ tvals, yvals ] = ode15s( @(t,y) sin(y)/t, [0 10], 1 );
[ tvals, xvals ] = ode15s( @(t,x) [ x(2); x(1)-4*x(2) ], [0 10], [1; 0] );

%----------------------------------------
% tic
% commands to time
% toc
disp( 'Timing for ode45' );
opt = odeset( 'RelTol', 1e-8 );
tic
[t,x] = ode45( @(t,x) [1000*x(2); -1000*x(1) - x(2)], [0 2000],...
               [-2;5], opt );
toc

%-------------------------------------------------------------------------------
%% Labs 5

%----------------------------------------
% plot( xvector, yvector )

% plot a phase portrait for a system given the solution xvals from
% ode45 in a new figure
figure();
plot( xvals(:,1), xvals(:,2), '-k', 'LineWidth', 2 );
% plot the component curves in a new figure
figure();
plot( tvals, xvals(:,1), '-k', tvals, xvals(:,2), '--k' );

%----------------------------------------
% plot3( xvec, yvec, zvec )
%     example command
plot3( xLS(:,1), xLS(:,2), xLS(:,3) );
title( 'Plot of Lorenz system trajectory' );
xlabel( 'x' ); ylabel( 'y' ); zlabel( 'z' );

%----------------------------------------
% PlotComet_3D( xvec, yvec, zvec, 'option', value,... );
%     example command
PlotComet_3D( xLS(:,1), xLS(:,2), xLS(:,3), 'Frequency', 50, 'blockSize', 100 );

%----------------------------------------
% copyright information
% (c) 2018, University of Michigan Mathematics Department
% licensed under a Creative Commons By-NC-SA 4.0 International License.
%----------------------------------------
