% Bare essentials of MATLAB for Math 46 class
% Alex Barnett 4/6/07, updated 4/3/08
%
% In windows or Mac, click on the Matlab icon.

% Anything after a % is a comment

% You can type each line at the command line, or save it as intro.m and
% run the whole thing as a program by typing `intro'

a = 5;                                % a symbol can be a number (scalar)
a                                     % without a ; the answer is printed
sa = sin(a)                           % sa now contains sin of a, ie sin(5)

b = [1 2 3]                           % a row vector (row vec)
c = [3; -2; 4.6]                      % a column vector (col vec)
c(3)                                  % look up the 3rd element of c
d = [1 2 3; 4 5 6]                    % 2x3 matrix (array)
d(2,1)                                % get entry at 2nd row, 1st col
d(:,1)                                % get whole 1st col
d(2,:)                                % get 2nd row
b*c                                   % matrix product. Careful with order!
% * computes matrix product. Therefore b*c computes dot product (a scalar),
% but c*b computes 3-by-3 `outer product'.
dot(b,c)                              % dot(c,b) same, better way for dot prod
d = c'                                % d is now transpose of c (row vec)
b.*d                                  % element-by-element multiply, a 1-by-3
sin(b)                                % sin of all elements at once (fast!)

plot(b, c, '+');                      % plot c vs b using + symbols
% Note: this means plot the 3 points (b(1),c(1)), (b(2),c(2)) and (b(3),c(3))
% b is interpreted as `x values' (horizontal), c as `y values' (vertical).
% Number of elements in b and c must match up (unless want multiple lines) 
axis([0 4 -3 5]);                     % manually choose axes range
plot(b, c, 'g-');                     % plot same using instead green lines

help                                  % list of help menus
help plot                             % give help on eg plot command
doc                                   % if you prefer graphical help window
whos                                  % show all currently-defined variables
save mydata                           % save all variables as file mydata.mat
quit                                  % need I say?
load mydata                           % load vars (can overwrite variables)
clear                                 % clear all variables (careful!)
% Use up and down arrows to scroll through recent commands, enter to execute
% I will postpone matrices & linear algebra for later...


% Solving general 1st-order ODE y' = f(t,y) numerically ----------------------
f = @(t,y) 1-y;                       % set up func of 2 vars f(t,y) = 1-y
f(2,3)                                % check it evaluates correctly

yo = 0;                               % initial condition at t=0 is yo=0
[ts, ys] = ode45(f, [0, 3], yo, odeset('reltol', 1e-10));      % solve ODE (using 'Runge-Kutta')
                                      % over the t-range [0,3]
% ts is now a col vec of t values, ys is col vec of corresponding y values
% ----------------------------------------------------------------------------

plot(ts, ys, '+-');                   % plot solution with + symbols and lines

yex = 1 - exp(-ts);                   % compute analytic (exact) y values
                                      % at each already-existing t value
% Note that we subtracted exp(ts) which is a col vec, from 1 which is a scalar,
% but it figured out to do the right thing (make 1 an appropriate vector of 1s)

hold on;                              % plot the next graph over this one
plot(ts, yex, 'ro');                  % compare exact solution as red o symbols
hold off;                             % stop overlaying future graphs

figure(2);                            % start a new figure
plot(ts, ys-yex, '+-');               % Subtract the two - best way to compare!
xlabel('t'); ylabel('y error');       % label axes: important to communicate!
title('dy/dt = 1-y: error in numerical solution ys-yex');


% Solving general 2nd-order ODE y'' = f(t,y,y') numerically ------------------
% eg consider u'' - u = eps.t.u
% rewrite as coupled 1st-order with col vec y = [u; v], where v = u'
% then ODE is y' = F(y) for F:R2->R2 some vector function

eps = 0.1;                               % note you have to set value of eps
F = @(t, y) [ y(2); y(1)*(1 + eps*t)];   % before you create vec func F
yo = [1;1];                              % IC for u and v
[ts, ys]  = ode45(F, [0 3], yo);         % numerically solve in t domain [0,3]
%-----------------------------------------------------------------------------

figure; plot(ts, ys(:,1), '-');          % plot just first col (u values) vs t
hold on; plot(ts, ys(:,2), '--r');       % add v = u' to the plot
xlabel('t');  legend('u', 'v = u'''); title('solution of 2nd-order ODE');
% note how in order to create a prime (') inside the text string I had to use
% '' (the third ' is the one which closes the string!)
% You could now compare against an analytic solution and add that to your plot


% See course website links for more complete introductions. Enjoy!