%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FIRST DAY MATLAB worksheet
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% for MTH 452-552
%%% this file helps to get familiar with basic elements of MATLAB
%%% for use in the Numerical Methods for ODE course
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%% M. Peszynska,  2005/1
%% Copyright Department of Mathematics, Oregon State University
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%%% INSTRUCTIONS: please type every line yourself
%%%               ignore all lines and text beginning with %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
diary mpesz
1 + 3
a = 1
b = 3
5 * a + b^2 + 1/a
1
1;
1,
%%% representation of numbers
format hex;1,format	%% use commas and semicolons
pi 
format long		%% to see more digits
pi
-2*pi+100+pi-100 	%% test if associative law applies
format short
-2*pi+100+pi-100 
eps			%% this is the machine accuracy
1+eps			%% that is, the smallest positive number added to 1 
			%% which makes the difference: 
			%% but you can see only that difference in ...
format hex
1 + eps
1 
1 + eps/2               %% eps/2 does not change the value
format
%%%%%%%%%%%%%%%%%%%%%% learn some programming: how to set up a loop
%%% compute factorial 
f=1; n=10; for i=1:n f=f*i;end;
f
f=1; n=150; for i=1:n f=f*n;end;f
f=1; n=200; for i=1:n f=f*i;end;f
%%% we have obtained an OVERFLOW: somewhere between 150! and 200!
format long
f
format hex;
f
format
%%%%%%%%%%%%%% EXERCISE: 
%%% find the smallest representable number: find n such that 1/n| is not zero

%%%%%%%%%%%%%%%%%%%%%% learn how to write a function saved in an M-file 
factorial(150)
%%%%%%%%%%%%%% EXERCISE: 
%%% write an M-file square.m that computes x^2: call it as square(pi)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% learn how to plot
y = 0:.1:1 %% note:	 	%this makes y into a vector
plot(y,sin(y))
y = 0:.01:10;
plot(y,sin(y));
plot(y,sin(y),'b',y,cos(y),'r');
y^2			%% y is a vector so you have to use different syntax
y.^2
%%%%%%%%%%%%%% EXERCISE: 
%%% plot various elementary functions in various intervals, e.g., square(y) ... 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% assess computational complexity: 
tic; f=1; n=150; for i=1:n f=f*i;end;f;toc
tic; factorial(150); toc
%%%%%%%%%%%%%% EXERCISE: 
%%% time the execution of 10000 times of computing pi^2 by calling your 
%%% function square.m and doing it directly

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% differential equations
%%% learn how to solve a differential equation symbolically
help dsolve
dsolve('Dy = -y + sin(t) ')
yfun = dsolve('Dy = -y + sin(t) ')
pretty(yfun)
yfun = dsolve('Dy = -y + sin(t) ','y(0)=1')
pretty(yfun)
%%%% how to plot direction fields
help quiver
[t,y] = meshgrid(0:.1:2,-1:.1:1);
u = ones(size(t)); v = funfile(t,y);
quiver (t, y, u, v)
%%% same thing with an inline function
funinline = inline('-y + sin(t)','t','y'); 
u = ones(size(t)); v = funinline(t,y);
quiver(t,y,u,v)
%%% numerical solution of ODEs: 
%%% there are many functions: here we use ode45
ode45(@funfile,[0,2],[1]);
ode45(funinline,[0,2],[1]);
ode45(@funfile,[0,2],[0]);	%% changed initial conditions
ode45(@funfile,[0,10],[1]);     %% plot in larger interval
tic;u = ones(size(t)); v = funinline(t,y);toc
tic;u = ones(size(t)); v = funfile(t,y);toc
%%% how to solve a system numerically
[t,y]=ode45(@harmonic,[0,2*pi],[1 0]);
plot(t,y(:,1))
plot(t,y(:,2))
plot(y(:,1),y(:,2))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% additional material
%%%%%%%%%% EXERCISE: play with fzero to solve x=exp(-x/2)
help fzero %%% try fzero to solve (numerically) the class example
%%%%%%%%%% EXERCISE: learn how to evaluate integrals numerically
help integral %%% no function of this name, try lookfor
lookfor integral 
help quad