Textbook: MATLAB, An Introduction with Applications, Amos Gilat (Fifth Edition)

Some Basics of MATLAB

Operations
clear
format short
%: use for description or comment
5+2 % Addition

ans = 7

7-3 % Subtraction

ans = 4

5*3 % Multiplication

ans = 15

7/5 % Right Division

ans = 1.4000

5\10 % Left Divison

ans = 2

4^3 % Exponentiation

ans = 64

3*5; % The output is not displayed

The order of operations:

1) Pharanthesis

2) Exponentiation

3) Multiplication and division

4) Addition and subtraction

5^1/3+7^0.3

ans = 3.4595

5^(1/3)+7^0.3

ans = 3.5028

Defining Scalar Variables and Some Math. Functions

x = pi % Create variables with the equal sign (=)

1
x = 3.1416

a=5*x+3

a = 18.7080

y = sin(-5) % You can provide inputs to a function using parentheses.

y = 0.9589

b=x+y

b = 4.1005

tan(x*y) % sin(x), cos(x), cot(x),tan(x). Here 'x' is in radians.

ans = -0.1298

z=-0.1298

z = -0.1298

c=(a-b)+40-x/y*z

c = 55.0327

abs(z) % Used for absolute value

ans = 0.1298

log(z) % Used for logarithm

ans = -2.0418 + 3.1416i

sqrt(x) % Used for square root

ans = 1.7725

nthroot(27,3) % Used for nth root

ans = 3

exp(x) % Used for exponential

ans = 23.1407

round(exp(x)) % Used for rounding to the nearest integer

ans = 23

cosd(60) % sind(x), cosd(x), cotd(x),tand(x). Here 'x' is in

[Link](4) % 4!

ans = 0.5000
ans = 24

Case1=z*b % You can also name a variable

2
 Case1 = -0.5322

ans % 'ans' gives the last answer that was not assgined to a specific
variable

ans = 24

inf % Used for infinity

ans = Inf

5+3i % i; imaginary part

ans = 5.0000 + 3.0000i

3-4j % j; imaginart part

ans = 3.0000 - 4.0000i

0/0 %NaN: Not a Number

ans = NaN

Managing Variables
x = pi

x = 3.1416

format long % Use for displaying 15 decimal digits

x

x =
3.141592653589793

format short % Use for displaying 4 decimal digits

x

x = 3.1416

format short g % Use for displaying 5 digits

12.6677846595998587

ans =
12.668

format long g % Use for displaying 15 digits

12.6677846595998587

ans =
12.6677846595999

format bank % Use for 2 decimal digits

12.6677846595998587

ans =
12.67

3
 290/7 % Follows the last format

ans =
41.43

digits(7) % Use for displaying any arbitrary digits. It is used with 'vpa'.
vpa(pi)

ans =

digits(15) % Use for displaying any arbitrary digits. It is used with 'vpa'.
vpa(290/7)

ans =

digits % Use for displaying the number of digits.

Digits = 15

who % displays the current variables in the memory

Your variables are:

Case1 a ans b c x y z

% clear a b % clears the variables a and b from the workspace

% clear % clears the all variables from the work space

Examples

clear x;
x=24

x = 24

tand(3*x)

ans = 3.0777

4
 (3*tand(x)-tand(x)^3)/(1-3*tand(x)^2)

ans = 3.0777

cosd(4*x)

ans = -0.1045

8*(cosd(x)^4-cosd(x)^2)+1

ans = -0.1045

R1=120.6;
R2=119.3;
R3=121.2;
R4=118.8;
V=14;
Vab=V*(R1*R3-R2*R4)/((R1+R2)*(R3+R4))

Vab = 0.1079

Homeworks: CH1-12, 26

pdffile = export("MATLAB CHAPTER [Link]","[Link]")

pdffile =
'/MATLAB Drive/[Link]'

5
Creating Arrays
One Dimensional Array (Vector)
clear
4 % scalar

ans = 4

r=[3 5] % row vector: Variable Name= [Vector Elements] without a

punctuation between them

r = 1×2
3 5

p=[1;3] % column vector: Variable Name= [Vector Elements] with semicolon (;)
between them

p = 2×1
1
3

g=[1:5] % variable name= [m:n] creates a vector from m to n spaced by 1

g = 1×5
1 2 3 4 5

h=[1:3:19] % variable name= [m:q:n] creates a vector from m to n spaced by q

h = 1×7
1 4 7 10 13 16 19

f=linspace(1,10,5) % variable name=linspace (m,q,n), where m: first element,

q:last element, n:number of elements. The vector is equally spaced.

f = 1×5
1.0000 3.2500 5.5000 7.7500 10.0000

Two Dimensional Array (Matrix)

A=[3 4 5;6 7 8;2 4 5] % matrix: Variable name=[1st row; 2nd row;...;last row]

A = 3×3
3 4 5
6 7 8
2 4 5

B=[linspace(10,30,6);(1:6);2 3 7 5 7 1]

B = 3×6
10 14 18 22 26 30
1 2 3 4 5 6
2 3 7 5 7 1

AA=2;
AB=3;
fe=4;

1
 C=[AA*fe AB+AA;fe*AA sqrt(fe)]

C = 2×2
8 5
8 2

rand(3) % rand(m) creates a square matrix with m rows and m columns

ans = 3×3
0.8147 0.9134 0.2785
0.9058 0.6324 0.5469
0.1270 0.0975 0.9575

zr=zeros(2,4) % zeros(m,n) creates a matrix with m rows and n columns and

all elements are 0

zr = 2×4
0 0 0 0
0 0 0 0

ne=ones(3,5) % zeros(m,n) creates a matrix with m rows and n columns and all
elements are 1

ne = 3×5
1 1 1 1 1
1 1 1 1 1
1 1 1 1 1

idn=eye(4) % eye(m) creates an identity matrix with m rows and m columns.

idn = 4×4
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1

M=[1 4;2 3];

Mtr=M' % quote(') creates a transpose of a matrix.

Mtr = 2×2
1 2
4 3

Array Addressing (Vector)

Vec1=[-5 5 13 2 6 24 9]; % Variable name=[row vector]
Vec1(2) % Variable name(n): chooses the nth element of the row vector.

ans = 5

Vec1(6)

ans = 24

Vec2=[2;3;5;7;21;12]; % Variable name=[column vector]

Vec2(4) % Variable name(n): chooses the nth element of the column vector.

ans = 7

2
 Vec2(6)

ans = 12

VV=Vec1(3)*sqrt(Vec2(2))

VV = 22.5167

Vec1

Vec1 = 1×7
-5 5 13 2 6 24 9

Vec1(2:5) % vn(m:n): chooses the elements m through n of the vector vn.

ans = 1×4
5 13 2 6

Vec1([3,5,7])

ans = 1×3
13 6 9

Vec1([3,5:7])

ans = 1×4
13 6 24 9

Array Addressing (Matrix)

Mat1=[(2:5);12 4 5 8; linspace(3,6,4);3 5 7 8] % Variable name=[matrix]

Mat1 = 4×4
2 3 4 5
12 4 5 8
3 4 5 6
3 5 7 8

Mat1(4,2) % Variable name(m,n): chooses the element in row m and column n.

ans = 5

Mat1(3,4)^Mat1(2,3)

ans = 7776

Mat1

Mat1 = 4×4
2 3 4 5
12 4 5 8
3 4 5 6
3 5 7 8

Mat1(:,2) % VN(:,m) refers to the entire mth column of the matrix VN

ans = 4×1
3
4

3
 4
5

Mat1(3,:) % VN(m,:) refers to the entire mth row of the matrix VN

ans = 1×4
3 4 5 6

Mat1(:,2:3) % VN(:,m:n) refers to the entire rows between columns m and n of

the matrix VN

ans = 4×2
3 4
4 5
4 5
5 7

Mat1(2:3,:) % VN(m:n,:) refers to the entire columns between rows m and n of

the matrix VN

ans = 2×4
12 4 5 8
3 4 5 6

Mat1(2:4,1:2) % VN(m:n:,p:q) refers to the elements in rows m through nand

columns p through q of the matrix VN.

ans = 3×2
12 4
3 4
3 5

Mat2=[ones(1,6);zeros(1,6);[0:3:15];linspace(4,9,6)]

Mat2 = 4×6
1 1 1 1 1 1
0 0 0 0 0 0
0 3 6 9 12 15
4 5 6 7 8 9

Mat3=Mat2([1,3],[1,5:6])

Mat3 = 2×3
1 1 1
0 12 15

Adding Elements to Existing Variables (to a Vector)

MX=[1:5]

MX = 1×5
1 2 3 4 5

MX(6)=9 % VN(n)=m assigns a value 'm' to the 'n'th element

MX = 1×6
1 2 3 4 5 9

4
 MX(5:10)=[10:5:35] % It adds 6 elements starting with the 5th.

MX = 1×10
1 2 3 4 10 15 20 25 30 35

C=[2 4 1]

C = 1×3
2 4 1

C(6)=5 % It assigns a value 5 to the 6th element. The others are assigned to
the zeros.

C = 1×6
2 4 1 0 0 5

S(7)=8 % It assigns a value 8 to the 7th element of a new vector. The

others are assigned to the zeros.

S = 1×7
0 0 0 0 0 0 8

R1=[linspace(4,16,4)];
R2=[1:4];
RR=[R1 R2] % If M is a row vector and N is another row vector, VN=[M,N] is a
new row vector appending M and N.

RR = 1×8
4 8 12 16 1 2 3 4

RRT=[R1';R2'] % VNT=[M',N'] generates a new column vector VNT by appending

M' and N'.

RRT = 8×1
4
8
12
16
1
2
3
4

Adding Elements to Existing Variables (to a Matrix)

D=[2 3 4; 2 5 8]

D = 2×3
2 3 4
2 5 8

D(3,:)=[15:4:23] % VN(M,:)=[Vec] adds the vector 'Vec' as the 'M'th row of

VN.

D = 3×3
2 3 4
2 5 8
15 19 23

5
 E=eye(3)

E = 3×3
1 0 0
0 1 0
0 0 1

F=[D E] % It appends the matrix E to matrix D. The numbers of rows in E and

D must be the same.

F = 3×6
2 3 4 1 0 0
2 5 8 0 1 0
15 19 23 0 0 1

D = 3×3
2 3 4
2 5 8
15 19 23

D(5,4)=12 % VN(m,n) assigns a value to the (m,n) element. The other elements
are assigned to the zeros.

D = 5×4
2 3 4 0
2 5 8 0
15 19 23 0
0 0 0 0
0 0 0 12

N(3,4)=45 % It assigns a value to the (3,4) element of a new matrix. The

other elements are assigned to the zeros.

N = 3×4
0 0 0 0
0 0 0 0
0 0 0 45

Deleting Elements from Existing Variables (from a Vector)

DD=[2:2:28]

DD = 1×14
2 4 6 8 10 12 14 16 18 20 22 24 26

DD(7)=[] % VN(m)=[] deletes the 'm'th element of the vector VN.

DD = 1×13
2 4 6 8 10 12 16 18 20 22 24 26 28

DD(3:6)=[] % VN(m:n)=[] deletes the elements m through n of the vector VN.

DD = 1×9
2 4 16 18 20 22 24 26 28

Deleting Elements from Existing Variables (from a Matrix)

6
 EE=[4 6 12 5 23; 4 7 23 1 56; 1:5]

EE = 3×5
4 6 12 5 23
4 7 23 1 56
1 2 3 4 5

EE(:,4:5)=[] % VN(:,m:n)=[] deletes all the columns m through n.

EE = 3×3
4 6 12
4 7 23
1 2 3

EE(:,2)=[] % VN(:,m)=[] deletes the column m.

EE = 3×2
4 12
4 23
1 3

EE(:,end)=[] % VN(:,end)=[] deletes the last column.

EE = 3×1
4
4
1

FF=[4 6 12 5 23; 4 7 23 1 56; 1:5; 34 5 7 78 9]

FF = 4×5
4 6 12 5 23
4 7 23 1 56
1 2 3 4 5
34 5 7 78 9

FF(2:4,:)=[] % VN(m:n,:)=[] deletes all the rows m through n.

FF = 1×5
4 6 12 5 23

GG=[4 6 12 5 23; 4 7 23 1 56; 1:5]

GG = 3×5
4 6 12 5 23
4 7 23 1 56
1 2 3 4 5

GG(2,:)=[] % VN(m,:)=[] deletes the row m.

GG = 2×5
4 6 12 5 23
1 2 3 4 5

Managing Arrays (Continued)

V=[2 4 1 5 6];
length(V) % length(M) gives the number of elements in the vector V

7
ans = 5

size(GG) % size(M) gives the size mxn of the array M and the result is a row
vector [m,n].

ans = 1×2
2 5

GG

GG = 2×5
4 6 12 5 23
1 2 3 4 5

reshape(GG,5,2) % reshape(M,m,n) creates an m by n matrix from the elements

of matrix A. The elements are taken column after column.

ans = 5×2
4 3
1 5
6 4
2 23
12 5

a=[1 4 6]; % diag(M) results in a diagonal matrix with the elements of M

(vector) in the diagonal.
ad=diag(a)

ad = 3×3
1 0 0
0 4 0
0 0 6

b=rand(4)

b = 4×4
0.9649 0.4854 0.9157 0.0357
0.1576 0.8003 0.7922 0.8491
0.9706 0.1419 0.9595 0.9340
0.9572 0.4218 0.6557 0.6787

bd=diag(b) % diag(M) gives a column vector from the diagonal elements of M

(matrix).

bd = 4×1
0.9649
0.8003
0.9595
0.6787

pdffile = export("MATLAB CHAPTER [Link]","[Link]")

pdffile =
'/MATLAB Drive/MATLAB CHAPTER 1/[Link]'

8
Mathematical Operations with Arrays
Addition
Va=[2,4,12,4];
Vb=[3,4,1,7];
Vc=Va+Vb % Adds the vectors Va and Vb. Vector sizes must agree.
Vc = 1×4
5 8 13 11
Vd=Va-Vb % Subtracts the vectors Va and Vb. Vector sizes must agree.
Vd = 1×4
-1 0 11 -3
Vc-3
ans = 1×4
2 5 10 8
MA=[1:4;linspace(2,8,4);4:3:13]
MA = 3×4
1 2 3 4
2 4 6 8
4 7 10 13
MB=[2,3,12,4;5:4:20;ones(1,4)]
MB = 3×4
2 3 12 4
5 9 13 17
1 1 1 1
MC=MA+MB % Adds the matrix MA and MB. Matrix sizes must agree.
MC = 3×4
3 5 15 8
7 13 19 25
5 8 11 14
MD=MA-MB % Subtracts the matrix MA and MB. Matrix sizes must agree.
MD = 3×4
-1 -1 -9 0
-3 -5 -7 -9
3 6 9 12
MC+4
ans = 3×4
7 9 19 12
11 17 23 29
9 12 15 18

Multiplication
M1=[2,3,5;1,4,6;2,12,9;4,6,3]
M1 = 4×3
2 3 5
1 4 6
2 12 9
4 6 3

1
 M2=[3,1;4,1;3,2]
M2 = 3×2
3 1
4 1
3 2
M3=M1*M2 % 'm x q * q x n = m x n': Inner matrix dimensions (q) must agree.
M3 = 4×2
33 15
37 17
81 32
45 16
V1=[2,4,5]
V1 = 1×3
2 4 5
V2=[3;2;6]
V2 = 3×1
3
2
6
V3=V1*V2 % Dot product of two vectors.
V3 = 44
a=4;
Ma=M1*a
Ma = 4×3
8 12 20
4 16 24
8 48 36
16 24 12

A system of linear equations (e.g. three equations; three unknowns):

Matrix form:

In matrix notation:

2
Inverse of a Matrix:
A=[3,5,6;2,4,5;-2,5,-4]
A = 3×3
3 5 6
2 4 5
-2 5 -4
B=inv(A) % inv(A) or A^(-1) is the inverse of the matrix A.
B = 3×3
1.6400 -2.0000 -0.0400
0.0800 0 0.1200
-0.7200 1.0000 -0.0800
A*B % A*A^(-1)=I : Identity Matrix
ans = 3×3
1.0000 0 0.0000
0 1.0000 0
0 0 1.0000

Determinants (e.g. 2x2 matrix):

det(A) % det(A) is the determinant of the matrix A.

ans = -25.0000

Division
Left division: is obtained by the left division:

in MATLAB.

Right division: is obtained by the right

division: in MATLAB.

Example: Solve the following system of linear equations:

3
 is equivalent to the matrix form , i.e.,

or , i.e.,

clear A B
A=[2,4,-3;5,-10,3;1,-1,2]
A = 3×3
2 4 -3
5 -10 3
1 -1 2
B=[3;2;4]
B = 3×1
3
2
4
X=A\B % Left division
X = 3×1
1.7273
1.1818
1.7273
x=inv(A)*B
x = 3×1
1.7273
1.1818
1.7273
C=[2,5,1;4,-10,-1;-3,3,2]
C = 3×3
2 5 1
4 -10 -1
-3 3 2
D=[3 2 4]
D = 1×3
3 2 4
clear X x
X=D/C % Right division
X = 1×3
1.7273 1.1818 1.7273
x=D*inv(C)
x = 1×3
1.7273 1.1818 1.7273

4
Array Operations (or Element-by-Element Operations)
Array operations (or element-by-element operations) can be done only with arrays of the same size.
clear A B C D
A=[1,2;3,2]
A = 2×2
1 2
3 2
B=[2,4;1,3]
B = 2×2
2 4
1 3
A.*B % The code .* is used for element-by-element multiplication.

A./B % The code ./ is used for element-by-element division.

ans = 2×2
0.5000 0.5000
3.0000 0.6667
A.^B % The code .^ is used for element-by-element exponentiation.
ans = 2×2
1 16
3 8
clear x y
x=[1:5]
x = 1×5
1 2 3 4 5
x.^2
ans = 1×5
1 4 9 16 25

Calculate:

y=(x.^3 + x)./(4.*x.^2)
y = 1×5
0.5000 0.6250 0.8333 1.0625 1.3000
sin(y) % gives the sinus of all the elements in array (vector) y
ans = 1×5
0.4794 0.5851 0.7402 0.8736 0.9636
sqrt(A) % gives the square root of all the elements in array (matrix) A
ans = 2×2
1.0000 1.4142
1.7321 1.4142

Some Array Functions

5
e=[2 10 3 4];
mean(e) % gives the mean value of the elements (e.g. (2+10+3+4)/4) of the
vector e.
ans = 4.7500
max(e) % gives the maximum value of the elements of the vector e.
ans = 10
min(e) % gives the minimum value of the elements of the vector e.
ans = 2
[d,n]=max(e) % d:the maximum value of e and n: the position of the maximum
value of e
d = 10
n = 2
[d,n]=min(e) % d:the minimum value of e and n: the position of the minimum
value of e
d = 2
n = 1
e
e = 1×4
2 10 3 4
sum(e) % gives the summation of all the elements of the vector e.
ans = 19
sort(e) % gives the ascending order of the elements of the vector e.
ans = 1×4
2 3 4 10
median(e) % gives the median value of the elements (e.g. (3+4)/2) of the
vector e.
ans = 3.5000
std(e) % gives the standard deviation of the elements of the vector e.
ans = 3.5940
clear a
a=[1,3,-2];
b=[2,-5,1];
dot(a,b) % dot(m,n) gives the dot products of two vectors (row or column) m
and n.
ans = -15
cross(a,b) % cross(m,n) gives the cross product (axb) of two vectors m and
n. The number of elements of each vector must agree.
ans = 1×3
-7 -5 -11
c=[2;-5;1];
dot(a,c)
ans = -15

6
 cross(a,c)
ans = 1×3
-7 -5 -11

Commands for Random Numbers

• rand Command:

rand % gives a random number btw 0 and 1.

ans = 0.4899
rand(1,5) % rand(1,n) gives n element row vector of random numbers btw 0 and
1.
ans = 1×5
0.6949 0.4114 0.0348 0.2928 0.8014
rand(4) % rand(n) gives n x n matrix with random numbers btw 0 and 1.
ans = 4×4
0.3465 0.7395 0.1895 0.0161
0.0833 0.5247 0.1237 0.8960
0.5111 0.8045 0.8210 0.5154
0.3668 0.8169 0.6379 0.5445
rand(3,4) % rand(m,n) gives m x n matrix with random numbers btw 0 and 1.
ans = 3×4
0.6064 0.3829 0.3320 0.8282
0.7604 0.0846 0.8397 0.1765
0.8553 0.7339 0.3717 0.1295
randperm(6) % randperm(n) gives n element row vector of integer random
numbers btw 1 and n.
ans = 1×6
2 6 5 3 4 1
randperm(5,3) % randperm(m,n) gives n element row vector of random numbers
btw 1 and m.
ans = 1×3
5 3 4
10*rand+5 % (b-a)*rand+a gives a random values btw a and b (e.g. btw a=5 and
b=15)
ans = 13.9450
h=12*rand(2,3)-2 %(b-a)*rand(m,n)+a gives a random values btw a and b (e.g.
btw a=-2 and b=10)
h = 2×3
0.9643 2.9064 -0.2763
1.7281 6.4961 8.4559
round(h)
ans = 2×3
1 3 0
2 6 8

7
 • randi Command:

randi(22) % randi(n) gives an integer random number btw 1 and n.

ans = 2
randi(25,4) % randi(m,n) gives an nxn matrix with integer random numbers btw
1 and m.
ans = 4×4
12 6 14 13
1 17 23 7
19 14 11 24
18 22 9 12
randi(32,3,4) % randi(a,m,n) gives an mxn matrix with integer random numbers
btw 1 and a.
ans = 3×4
9 13 19 8
14 6 12 17
23 28 8 14
randi([12,45],2,4) % randi([a,b],m,n) gives an mxn matrix with integer
random numbers btw a and b.
ans = 2×4
37 40 16 18
14 35 41 32

• randn Command

randn(3,4) % rand(m,n) gives an mxn matrix with integer random numbers with
mean 0 and standard deviation of 1.
ans = 3×4
0.2195 0.5356 0.4189 0.1509
-0.8782 -0.3420 -0.6897 -0.9913
-2.8614 -0.5979 0.6712 0.8368
l=3*randn(1,5)+20
l = 1×5
21.4235 23.7565 17.3221 17.3142 20.9400
round(l)
ans = 1×5
21 24 17 17 21

8
PLOTS
% Plot of a given data:
r=[0:5];
d=[3 7 12 20 24 33];
plot(r,d) % plot(x,y) generates two-dimensional plots. Here x and y are
vectors.

Line Specifiers
• Line style: - (solid) -- (dashed) : (dotted) -. (dash-dot)

plot(r,d,'--') % plot(x,y,'line style')

1
 • Line color: r (red) , g (green) , b (blue) , c (cyan) , m (magenta) , y (yellow) , k (black) , w (white)

plot(r,d,':r') % plot(r,d,'line style line color')

2
 • Market type: + , o , * , . , x , ^ , v , s , d , p , h , < , >

plot(r,d,'-.kh') % plot(r,d,'line style line color market type') The order

of the specifiers is not important.

3
plot(r,d,'r-.*') % plot(r,d,'line style line color market type') The order
of the specifiers is not important.

4
plot(r,d,'-gv') % plot(r,d,'line style line color market type')

5
plot(r,d,'b') % plot(r,d,'line style line color market type')

6
plot(r,d,'b*') % plot(r,d,'line style line color market type')

7
Property Name and Property Value
• 'LineWidth' determines the width of the line.
• 'MarkerSize' determines the size of the marker.
• 'MarkerEdgeColor' determines the color of the edge of the marker.
• 'MarkerFaceColor' determines the color of the filling for filled markers.

plot(r,d,'r-.*','LineWidth',2,'MarkerSize',10,'MarkerEdgeColor','k') %
plot(r,d,'LineSpecifiers','PropertyName',PropertyValue)

8
plot(r,d,'r:>','LineWidth',2,'MarkerSize',10,'MarkerEdgeColor','g','MarkerFac
eColor','y')

9
% Plot of a function
x=linspace(1,20,8);
y=sin(3*x).*cos(2*x)+x.^2;
plot(x,y,'--ko','MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',10)

10
% Plot of a function:
fplot('5*(t^2)+10*t+5+cos(t)',[2 15],'-ko','LineWidth',1) % fplot('function',
[xmin xmax],'Line Specifiers')
Warning: fplot will not accept character vector or string inputs in a future release. Use
fplot(@(t)5.*(t.^2)+10.*t+5+cos(t)) instead.

11
Plotting Multiple Graphs in the Same Plot
• By the plot Command

clear x y
x=[-8:0.01:2];
y=3*x.^2+6*x+exp(2*x)+5;
yd=6*x+6+2*exp(2*x); % First derivative of y
ydd=6+4*exp(2*x); % Second derivative of y
plot(x,y,'--r',x,yd,':b',x,ydd,'-.k','LineWidth',2) % plot(first plot,second
plot,third plot,etc...,options)

12
 • By the hold on and hold off commands

plot(x,y,'--r','LineWidth',2)
hold on
plot(x,yd,':b','LineWidth',4)
plot(x,ydd,'-.k','LineWidth',1)
hold off

13
yddd=8*exp(2*x);
plot(x,yddd,'-g','LineWidth',2)

14
 • By the line command

h=[-2:2];
g=5*h.^2+5*h;
plot(h,g,'-ko','LineWidth',2,'MarkerSize',10,'MarkerFaceColor','y')
line(x,y,'LineStyle','-','Color','r','LineWidth',2) %
line(x,y,'PropertyName',PropertyValue) does not have 'LineSpecifiers'
line(x,yd,'LineStyle',':','Color','b','LineWidth',2)

15
Formatting Plot

16
%Formatting Plot by Commands
plot(x,y,'--r',x,yd,':b',x,ydd,'-.k','LineWidth',2)
xlabel('\it{Interval of x}','Color','r') % xlabel('The name of the x-
axis',options)
ylabel('\it{y(x) and its derivatives}','Color','g') % ylabel('The name of
the y-axis',options)
title('\bf{y,yd,ydd vs x}','FontSize',20) % title('The name of the
plot',options)
%legend('y vs x','yd vs x','ydd vs x','LineWidth',3,'EdgeColor','c') %
legends('The name graph 1','The name graph 2','others', options)
%legend('y vs x','yd vs x','ydd vs
x','LineWidth',3,'EdgeColor','c','Location','northwest')

17
 legend('y^{2\alpha} vs x_{\beta}','yd vs x','ydd vs
x','LineWidth',3,'EdgeColor','c','Location','northwest','NumColumns',2)
text(-2,100,'x_{1}=Coord.1','FontSize',15,'Color','m') % text(x,y,'name of
that specific coord.')
text(1,150,'x_{2}=Coord.2','FontSize',15,'Color','m')
% axis equal : Both axes are set to the same scale
% axis square : The area around the axes is squared.
% axis tight: The axes limits are ordered for the data range.
axis([-5,2,-25,160]) % axis([xmin,xmax,ymin,ymax]) changes the limits of the
axes x and y.
grid on % Adds the grid lines to the plot. 'grid off' command removes the
grid lines from the plot.

Error Bars
• Symmetric error bar:

x=linspace(1,20,8);
y=sin(3*x).*cos(2*x)+x.^2;
plot(x,y,'--ko','MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',5)
yer=[5,12,7,23,21,12,24,18]; % errors at each point of the data
hold on

18
errorbar(x,y,yer) % errorbar(x,y,e), e is a vector consisting of the value
of the errors. The sizes of the vectors x,y,and e must agree.
yerh=[1,1,0.5,0.3,0.8,1,2,0.2];
errorbar(x,y,yerh,'horizontal')
hold off

• Nonsymmetric error bar:

x=linspace(1,20,8);
y=sin(3*x).*cos(2*x)+x.^2;
plot(x,y,'--ko','MarkerEdgeColor','k','MarkerFaceColor','g','MarkerSize',5)
yu=[5,12,11,7,15,9,9,8]; % upper limits of the erros at each point
yl=[4,7,10,12,9,13,14,5]; % lower limits of the erros at each point
hold on
errorbar(x,y,yu,yl) % errorbar(x,y,u,l) u is a vector consisting of the
upper limits of the errors and l is a vector consisting of the lower limits
of the errors. The sizes of the vectors x,y,u,and l must agree.

19
Special Graphs
• Vertical Bar Plots

yr=[2011:2019];
NoPO=[20 30 10 25 20 45 40 23 30];
bar(yr,NoPO,'g','BarWidth',0.4) % bar(x,y,options) gives vertical bar plots.
xlabel('year')
ylabel('The Number of Police Officers')
title('The Number of Police Officers in Crimeville, 2011 to 2019')

20
 • Horizontal Bar Plots

yr=[2011:2019];
NoPO=[20 30 10 25 20 45 40 23 30];
barh(yr,NoPO,'g','BarWidth',0.4) % barh(x,y,options) gives horizontal bar
plots.

21
 • Stairs Plot

yr=[2011:2019];
NoPO=[20 30 10 25 20 45 40 23 30];
stairs(yr,NoPO,'b') % stairs(x,y,options) gives stairs plots.

22
 • Stem Plot

yr=[2011:2019];
NoPO=[20 30 10 25 20 45 40 23 30];
stem(yr,NoPO,'b') % stem(x,y,options) gives stem plots.

23
 • Pie Plot

NoPO=[5 30 20 45 40 10];
Cities = {'A','B','C','D','E','F'};
pie(NoPO,Cities) % pie(x) gives the relative size of different values. %
pie(x,label) assignes labels to the x vector.
title('NoPO in Cities')

24
HISTOGRAMS
Plots that display the distribution of data are called histograms. The histogram displays the number of data
points in each subranges, or bins. The whole range of a particular collection of data points is split into bins.
Gr=[23 25 12 67 54 89 83 76 34 92 34 89 42 67 98 18 81 32 67 89 77 66 76 69
48 89 91 78]; % The grades of the 28 students
hist(Gr) % the number of data points in each bin is plotted with the range
of data divided into 10 equally spaced bins.
xlabel('\bf{Grades}')
ylabel('\bf{Number of Students}')
title('Histogram of Grades Data','FontSize',15,'Color','r')

25
The smallest value in 'Gr' is 12 and the largest values in 'Gr' is 98. The range of the data is the 86 (98-12). The
range of the data is divided into 10 equally spaced subranges. The width of each subrange is then 8.6 (86/10).
Therefore, the range of the first subrange is from 12 to 20.6 (12+8.6) and contains two points (12 and 18) ; the
range of the second subrange is from 20.6 to 29.2 (20.6+8.6) contains two points (23 and 25); and so on. One
of the subrange from 55 to 63.6 does not contain any points.
Gr=[23 25 12 67 54 89 83 76 34 92 34 89 42 67 98 18 81 32 67 89 77 66 76 69
48 89 91 78]; % The grades of the 28 students
hist(Gr,6) % hist(y,nbins) divides the range in 'n' equally spaced subranges.

26
The smallest value in 'Gr' is 12 and the largest values in 'Gr' is 98. The range of the data is the 86 (98-12). The
range of the data is divided into 6 equally spaced subranges. The width of each subrange is then 14.3 (86/6).
Therefore, the range of the first subrange is from 12 to 26.3 (12+14.3) and contains four points (12 ,18, 23, and
25) ; the range of the second subrange is from 26.3 to 40.6 (26.3+14.3) contains three points (34, 34, and 32);
and so on. One of the subrange from 55 to 63.6 does not contain any points.
Gr=[23 25 12 67 54 89 83 76 34 92 34 89 42 67 98 18 81 32 67 89 77 66 76 69
48 89 91 78]; % The grades of the 28 students
x=[20:15:95];
hist(Gr,x) % divides the range so that the center of each subrange is
locacated at the values of the x vector.

27
s10=hist(Gr) % gives the number of students in each subrange. The range is
divided into 10 equally spaced.
s10 = 1×10
2 2 3 1 2 0 5 4 6 3
s6=hist(Gr,6)% gives the number of students in each subrange. The range is
divided into 6 equally spaced.
s6 = 1×6
4 3 3 5 6 7
sx=hist(Gr,x) % gives the number of students in each subrange centered at
each value of x vector.
sx = 1×6
4 4 2 5 6 7
[s10,s10out]=hist(Gr) % s10out displays the value of the center of each
subrange.
s10 = 1×10
2 2 3 1 2 0 5 4 6 3
s10out = 1×10
16.3000 24.9000 33.5000 42.1000 50.7000 59.3000 67.9000 76.5000
[s6,s6out]=hist(Gr,6) % s6out displays the value of the center of each
subrange.
s6 = 1×6

28
 4 3 3 5 6 7
s6out = 1×6
19.1667 33.5000 47.8333 62.1667 76.5000 90.8333
[sx,scout]=hist(Gr,x) % scout displays the value of the center of each
subrange.
sx = 1×6
4 4 2 5 6 7
scout = 1×6
20 35 50 65 80 95

Polar Plots
t=[0:0.01:pi];
f=3*cos(0.5*t).^2+t;
polar(t,f) % polar(theta,radius(theta),options). theta and radius are
vectors.

Multiple plots on the same page

subplot command gives the plots on the same page.

subplot(m,n,p) creates m x n ordered figure page and p determines the location of a plot in this figure page.
x=linspace(1,20,8);

29
y=sin(3*x).*cos(2*x)+x.^2;
subplot(3,1,1)
plot(x,y)
h=[-2:2];
g=5*h.^2+5*h;
subplot(3,1,2)
plot(h,g)
r=[0:5];
d=[3 7 12 20 24 33];
subplot(3,1,3)
plot(r,d)

30
POLYNOMIALS
Values, Roots, and Coefficients of a Polynomial
Polynomial of degree n:

Matlab representation: p = [ ]

The value of a polynomial at x: polyval (p , x) or polyval ([ ] , x)

The roots of a polynomial: roots (p) or roots ([ ])

Determining the coefficients of a polynomial for given roots: poly(r) , where r is a vector of roots

Example

Consider a polynomial

a) Calculate

b) Plot in the domain

c) Find the roots of the polynomial

d) Confirm the coefficients in the polynomial by using poly command

Solution a)
format bank
p=[0.1 -0.2 -1 0 -41.5 235];
polyval(p,3) % or polyval([0.1 -0.2 -1 0 -41.5 235],3)
ans =
91.60

Solution b)
x=[-6:0.1:6];
y=polyval(p,x);
plot(x,y)

1
Solution c)
r=roots(p)
r = 5×1 complex
-
-
-

Solution d)
c=poly(r)
c = 1×6
1.00 -2.00 -10.00 0.00 -415.00

Example

a) Find the roots of a polynomial .

b) Confirm the coefficients by using poly command.

Solution a)
% p=[2 3 -11 -9 5]
% polyval(p) or

2
 r=roots([2 3 -11 -9 5])
r = 4×1
-2.78
2.02
-1.13
0.39

Therefore, the same polynomial may be rewritten in the following form:

Solution b)
c=poly(r)
c = 1×5
1.00 1.50 -5.50 -4.50 2.50

Operations with Polynomials

Addition:
format bank
p1=[0.1 -0.2 -1 0 -41.5 235]; % a polynomial of degree 5
p2=[3 -11 -9 5]; % a polynomial of degree 3
p=p1+[0 0 p2] % the shorter vector has to be the same size as the longer
vector. Therefore, two zeros are added to the shorter vector.
p = 1×6
0.10 -0.20 2.00 -11.00 -50.50

which implies the following answer:

Multiplication:
p12=conv(p1,p2) % The sizes of the polynomials do not have to be the same.
p12 = 1×9
0.30 -1.70 -1.70 13.30 -116.50

which implies the following answer:

Division:

Use the command [a , b] = deconv (p , q) for the division of the polynomials

p is the coefficients of the numerator polynomial

q is the coefficients of the denominator polynomial

a is the coefficients of the quotient polynomial

b is the coefficients of the remainder polynomial

Example

3
Divide the polynomial by .

p3=[2 -13 0 75 2 0 -60];

p4=[1 0 -5];
[a,b]=deconv(p3,p4)
a = 1×5
2.00 -13.00 10.00 10.00 52.00
b = 1×7
0 0 0 0 0

which implies the following answer: .

Derivatives:

The derivative of a polynomial p: polyder (p)

The derivative of a product of two polynomials p1 and p2: polyder (p1 , p2)

The derivative of a quotient of two polynomials p1 and p2: [a , b] = polyder (p1 , p2)

Example

Consider two polynomials by .

a) Calculate the derivative of and

b) Calculate the derivative of product *

c) Calculate the derivative of division

Solution a)
P1=[2 -6 3 1 -7];
P2=[1 0 -5];
d1=polyder(P1)
d1 = 1×4
8.00 -18.00 6.00 1.00

Answer:

d2=polyder(P2)
d2 = 1×2
2.00 0

Answer:

4
Solution b)
d3=polyder(P1,P2)
d3 = 1×6
12.00 -30.00 -28.00 93.00 -44.00

Answer:

Solution c)
[a,b]=polyder(P1,P2)
a = 1×6
4.00 -6.00 -40.00 89.00 -16.00
b = 1×5
1.00 0 -10.00 0 25.00

Answer: .

CURVE FITTING (with Polynomials)

The command: p = polyfit(x , y , n)

x: a set of points in the x-axis (independent variable)

y: a set of points in the y-axis (dependent variable)

n: the degree of the polynomial that fits the data

p: the coefficients of the polynomial that fits the data

Example:

Consider a set of seven points: (0.9, 0.9) ; (1.5, 1.5) ; (3, 2.5) ; (4, 5.1) ; (6, 4.5) ; (8, 4.9) ; (9.5, 6.3)
format bank
x=[0.9 1.5 3 4 6 8 9.5];
y=[0.9 1.5 2.5 5.1 4.5 4.9 6.3];
plot(x,y,'ro')

5
a) Find the coefficients of the polynomial for n=1 that fits the data and plot this polynomial.
p1=polyfit(x,y,1) % the degree of the polynomial is 1
p1 = 1×2
0.57 1.00

Answer:

xp=[0.9:0.1:9.5];
yp1=polyval(p1,xp); % the values of the polynomial p1 at each point xp
plot(x,y,'o',xp,yp1)

6
b) Find the coefficients of the polynomial for n=2 that fits the data and plot this polynomial.
p2=polyfit(x,y,2) % the degree of the polynomial is 2
p2 = 1×3
-0.06 1.20 -0.06

Answer:

xp=[0.9:0.1:9.5];
yp2=polyval(p2,xp); % the values of the polynomial p2 at each point xp
plot(x,y,'o',xp,yp2)

7
c) Find the coefficients of the polynomial for n=3 that fits the data and plot this polynomial.
p3=polyfit(x,y,3) % the degree of the polynomial is 3
p3 = 1×4
0.02 -0.40 2.61 -1.42

Answer:

xp=[0.9:0.1:9.5];
yp3=polyval(p3,xp); % the values of the polynomial p3 at each point xp
plot(x,y,'o',xp,yp3)

8
d) Find the coefficients of the polynomial for n=4 that fits the data and plot this polynomial.
p4=polyfit(x,y,4) % the degree of the polynomial is 4
p4 = 1×5
0.01 -0.19 1.06 -1.07 1.16
xp=[0.9:0.1:9.5];
yp4=polyval(p4,xp); % the values of the polynomial p4 at each point xp
plot(x,y,'o',xp,yp4)

9
e) Find the coefficients of the polynomial for n=5 that fits the data and plot this polynomial.
p5=polyfit(x,y,5) % the degree of the polynomial is 5
p5 = 1×6
-0.00 0.08 -0.84 3.78 -5.94
xp=[0.9:0.1:9.5];
yp5=polyval(p5,xp); % the values of the polynomial p5 at each point xp
plot(x,y,'o',xp,yp5)

10
f) Find the coefficients of the polynomial for n=6 that fits the data and plot this polynomial.
p6=polyfit(x,y,6) % the degree of the polynomial is 6
p6 = 1×7
-0.01 0.16 -1.79 9.39 -23.95

Answer:

xp=[0.9:0.1:9.5];
yp6=polyval(p6,xp); % the values of the polynomial p6 at each point xp
plot(x,y,'o',xp,yp6)

11
CURVE FITTING (with Functions other than Polynomials)
Other than polynomials, you can also fit the following functions to a set of given data.

All these functions may be written in the form of which can be fitted with a linear polynomial of
degree 1 (i.e. n=1).

12
Example:

Solution:
t=[0:0.5:5];
w=[6.00 4.83 3.70 3.15 2.41 1.83 1.49 1.21 0.96 0.73 0.64];
plot(t,w,'ko')

13
 •
Cannot be a power function ( ) since when .
•
Cannot be a logarithmic function ( ) since is in the data.
•
Can be an exponential function
•
Can be a reciprocal function ( )

Checking which one gives a straight line data points.

subplot(2,1,1)
plot(t,log(w),'ko') % w=b*exp(m*t) log(w)=(m*t)+b (line equation)
xlabel('t')
ylabel('log(w)')
title('in the form of log(w)=mt+b')
subplot(2,1,2)
plot(t,1./w,'ko') % w=1/((m*t)+b) (1/w)=(m*t)+b (line equation)
xlabel('t')
ylabel('1/w')
title('not in the form of w^{-1}=mt+b')

14
 •
Cannot be a reciprocal function ( )

•
Can be an exponential function

format bank
pe=polyfit(t,log(w),1)
pe = 1×2
-0.46 1.79

Answer:

tm=0:0.1:5;
wm=exp(1.79)*exp(-0.46*tm);
plot(t,w,'o',tm,wm)

15
Interpolation
Used for the prediction of values from a given set of data.

Curve fitting and interpolation are not the same thing. We create a curve through the data points during
interpolation. By doing this, we implicitly assume that the data points are precise and distinct. Curve fitting is
used when data include dispersion (noise), typically as a result of measurement errors.

• Linear Interpolation: yq = interp1(x , y , xq , ' linear ') or yq = linear(x , y , xq)

x = 0:pi/4:2*pi;
y = cos(x);
xq = 0:pi/16:2*pi; % query (or interpolation) points (You want to predict
the y values at these points)
yq1 = interp1(x,y,xq,'linear'); % interpolated values at the query points
plot(x,y,'o',xq,yq1,':.','MarkerSize',15)

16
 • Spline Interpolation: yq = interp1(x , y , xq , ' spline ') or yq = spline(x , y , xq)

Spline interpolation gives low-degree polynomials between each subsequent data points instead of giving a
single, high-degree polynomial to all of the values at once, i. e., '' the interpolant is a special type of piecewise
polynomial''. For instance (see the plot of the following example), instead of fitting a single polynomial of degree
eight to all the pairings of nine points, fit eight cubic polynomials between each pair.
yq2 = interp1(x,y,xq,'spline');
plot(x,y,'o',xq,yq2,':.','MarkerSize',15)

17
 • Pchip Interpolation: yq = interp1(x , y , xq , ' pchip ') or yq = pchip(x , y , xq)

x = -3:3;
y = [-1 -1 -1 0 1 1 1];
t = -3:0.01:3;
p = interp1(x,y,t,'pchip'); % or pchip(x,y,t)
s = interp1(x,y,t,'spline'); % or spline(x,y,t)
plot(x,y,'o',t,p,'r-',t,s,'k-.','MarkerSize',10,'LineWidth',2)
legend('data','pchip','spline')

18
pchip (Shape-Preserving Piecewise Cubic Interpolation) interpolates using a piecewise cubic polynomial
P(x) with these properties:

•
On each subinterval , the polynomial P(x) is a cubic Hermite interpolating polynomial
for the given data points with specified derivatives (slopes) at the interpolation points.
•
P(x) interpolates y, that is, , and the first derivative is continuous. The second

derivative is probably not continuous so jumps at the are possible.

•
The cubic interpolant P(x) is shape preserving. The slopes at the are chosen in such a way that P(x)
preserves the shape of the data and respects monotonicity. Therefore, on intervals where the data is
monotonic, so is P(x), and at points where the data has a local extremum, so does P(x).

spline constructs S(x) in almost the same way pchip constructs P(x). However, spline chooses the slopes
at the data points differently, namely to make even S˜(x) continuous. This difference has several effects:

• spline produces a smoother result, such that S˜(x) is continuous.

• spline produces a more accurate result if the data consists of values of a smooth function.

19
• pchip has no overshoots and less oscillation if the data is not smooth.
• pchip is less expensive to set up.
• The two are equally expensive to evaluate.

20
ORDINARY DIFFERENTIAL EQUATIONS
Example 1 (on First Order ODE): The velocity v of a probe falling towards Mars is governed by the differential
equation

where , , v is in , and t is in seconds. The probe enters the

atmosphere with an initial velocity of .

a) Find and plot the probe's velocity from to seconds.

tspan = [0 120]; % time interval

v0 = 6000; % initial velocity
[t,v] = ode45(@fallingbody,tspan,v0) % ODE function is defined at the bottom
of the script.
t = 57×1
0
0.2545
0.5090
0.7636
1.0181
1.7041
2.3900
3.0760
3.7620
4.6385

v = 57×1
103 ×
6.0000
5.7130
5.4523
5.2145
4.9966
4.4902
4.0784
3.7381
3.4503
3.1398

plot(t,v)
% Label axes
xlabel("time (seconds)")
ylabel("velocity (m/s)")

1
b) Find the probe's velocity after 120 seconds. Assign it to the variable v120.

v120=v(end)
v120 = 380.6347

Example 2 (on First Order ODE): The temperature of a pie cooling to room temperature after being taken out
of an oven is given by

where . When the pie is removed from the oven ( minute), the pie's temperature is

a) Find and plot the pie's temperature from to seconds.

tspan = [0 100]; % time interval

T0 = 175; % initial temperature in the oven
[t,T] = ode45(@pieODEfun,tspan,T0) % ODE function is defined at the bottom
of the script.

2
t = 45×1
0
2.0657
4.1314
6.1971
8.2628
10.7628
13.2628
15.7628
18.2628
20.7628

T = 45×1
175.0000
166.4578
158.3956
150.7866
143.6052
135.4515
127.8490
120.7606
114.1515
107.9890

plot(t,T)
% Label axes
xlabel("time (minutes)")
ylabel("temperature (\circC)")

3
b) Find and plot the pie's temperature from to seconds when the pie was removed from the

refrigerator ( minutes) instead of the oven. The pie's temperature is initially .

tspan = [0 100]; % time interval

T0 = 5; % initial temperature in the refregirator
[t,T] = ode45(@pieODEfun,tspan,T0); % ODE function is defined at the bottom
of the script.
plot(t,T)
% Label axes
xlabel("time (minutes)")
ylabel("temperature (\circC)")

4
Example 3 (on System of ODEs): Consider the following system of ODEs

The initial conditions are at .

a) Find and plot from to minutes.

tspan = [0 3]; % time interval

Y0 = [100,0,0];
[t,Y] = ode45(@threeODEs,tspan,Y0) % ODE function is defined at the bottom
of the script.
t = 81×1
0
0.0000
0.0000

5
 0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000

Y = 81×3
100.0000 0 0
99.9999 0.0001 0.0000
99.9999 0.0001 0.0000
99.9998 0.0002 0.0000
99.9998 0.0002 0.0000
99.9995 0.0005 0.0000
99.9993 0.0007 0.0000
99.9990 0.0010 0.0000
99.9988 0.0012 0.0000
99.9975 0.0025 0.0000

plot(t,Y(:,1),"r")
hold on
plot(t,Y(:,2),"b")
plot(t,Y(:,3),"m")
hold off
% Annotate plot
legend("x","y","z")
xlabel("time")
ylabel("Solutions to x, y, z")

6
b) Extract the column vectors x, y, z from the solution matrix Y. Assign them to the variables x, y, z,
respectively.

x=Y(:,1)
x = 81×1
100.0000
99.9999
99.9999
99.9998
99.9998
99.9995
99.9993
99.9990
99.9988
99.9975

x(end) % the last value of x

ans = 0.2479
min(x) % the minimum value of x
ans = 0.2479
y=Y(:,2)
y = 81×1
0

7
 0.0001
0.0001
0.0002
0.0002
0.0005
0.0007
0.0010
0.0012
0.0025

z=Y(:,3)
z = 81×1
0
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000

c) What's the maximum amount of y, and at what time does this maximum occur?
[ymax,imax] = max(y)
ymax = 51.9052
imax = 49
tmax = t(imax)
tmax = 0.6819

Example 4 (on Higher Order ODE): Consider the damped oscillatory motion, which is described by the
following second-order ODE:

where , , and . The spring is initially stretched

away from its unstretched length and the block is released from rest at that position.

a) Find the solution matrix (where the first and second columns are the solutions to the positions and velocities

of the mass, respectively) in the time from to in seconds.

The second-order ODE can be written in the following system of first-order ODEs.

8
 tspan = [0 1]; % time interval in seconds
Y0 = [0.1;0]; % initial conditions: y0=0.1 v0=0
[t,Y] = ode45(@harmonicMotion,tspan,Y0)
t = 189×1
0
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000

Y = 189×2
0.1000 0
0.1000 -0.0001
0.1000 -0.0001
0.1000 -0.0002
0.1000 -0.0002
0.1000 -0.0005
0.1000 -0.0007
0.1000 -0.0010
0.1000 -0.0012
0.1000 -0.0025

b) Extract the position of the mass from the first column of Y and store it in a variable y then plot y versus t.

y=Y(:,1);
plot(t,y)
% Label axes
xlabel("Time")
ylabel("Position")

9
c) Extract the velocity of the mass from the second column of Y and store it in a variable v then plot v versus t.

v=Y(:,2);
plot(t,v)
% Label axes
xlabel("Time")
ylabel("Velocity")

10
ODE function of Example 1:
function dvdt = fallingbody(t,v)
D = 3.3e-5;
g = 3.72;
dvdt = g - D*v^2;
end

ODE function of Example 2:

function dTdt = pieODEfun(t,T)
c = -0.028;
dTdt = c*(T-23);
end

ODE function of Example 3:

function dYdt = threeODEs(t,Y)
x = Y(1); % x
y = Y(2); % y
z = Y(3); % z
dxdt = -2*x;
dydt = 2*x - 0.02*y^2;
dzdt = 0.01*y^2;
dYdt = [dxdt;dydt;dzdt];

11
 end

ODE function for Example 4:

function dYdt = harmonicMotion(t,Y)

%Define constants k, m, and b

k = 50;
m = 0.05;
b = 0.5;

% TODO - Extract the position y from the first element of Y.

y = Y(1);

% TODO - Extract the velocity v from the second element of Y.

v = Y(2);

dydt = v;
dvdt = -(k/m)*y - (b/m)*v;

% TODO - Create dYdt, column vector containing dydt, and dvdt

dYdt = [dydt;dvdt];

end

12
SOLVING AN EQUATION WITH ONE VARIABLE
Example 1: Determine the solution of the equation .

format bank
fplot('exp(0.3*x)-x^2+4',[-20 20])
Warning: fplot will not accept character vector or string inputs in a future release. Use
fplot(@(x)exp(0.3.*x)-x.^2+4) instead.

x1=fzero('exp(0.3*x)-x^2+4',-2) % find the zero of function near -2

x1 =
-2.13

Answer: The value of the function is 0 at x=-2.13

x2=fzero('exp(0.3*x)-x^2+4',20) % find the zero of function near 20
x2 =
19.91

Answer: The value of the function is 0 at x=19.91

[x1,f1]=fzero('exp(0.3*x)-x^2+4',-2) % find the zero of function near -2
x1 =
-2.13
f1 =
0

1
 [x2,f2]=fzero('exp(0.3*x)-x^2+4',20) % find the zero of function near 20
x2 =
19.91
f2 =
-0.00

Example 2: Determine the solution of the equation .

fplot('cos(x)^2-3*x',[-2 2])
Warning: fplot will not accept character vector or string inputs in a future release. Use
fplot(@(x)cos(x).^2-3.*x) instead.

x=fzero('cos(x)^2-3*x',0.5) % find the zero of function near 0.5

x =
0.30

Answer: The value of the function is 0 at x=0.3

[x,f]=fzero('cos(x)^2-3*x',0.5) % find the zero of function near 0.5
x =
0.30
f =
0

FINDING A MINIMUM OR A MAXIMUM OF A FUNCTION

format bank

2
 fplot('exp(0.3*x)-x^2+4',[-20 20])
Warning: fplot will not accept character vector or string inputs in a future release. Use
fplot(@(x)exp(0.3.*x)-x.^2+4) instead.

[xmin,f]=fminbnd('exp(0.3*x)-x^2+4',10,20)
xmin =
15.45
f =
-131.67

Answer: The minimum value of the function is -131.67 at x=15.45

[xmax,f]=fminbnd('-exp(0.3*x)+x^2-4',-6,6) % multiply the function with
minus sign
xmax =
0.16
f =
-5.02

Answer: The maximum value of the function is 5.02 at x=0.16.

NUMERICAL INTEGRATION
Command: quad(function, a, b) or quadl(function, a, b)

Example: Use numerical integration to calculate the following integrals:

3
quad('cos(x).^2-3*x',-3,2)
ans =
9.74
quad('x.*exp(0.3*x)-x.^2+4',0,7)
ans =
24.59

4
THREE DIMENSIONAL PLOTS
Line Plots

The number of the elements of x, y, and z must be the same.

t=0:0.1:6*pi;
x1=sqrt(t).*sin(2*t);
y1=sqrt(t).*cos(2*t);
z1=0.5*t; % x1(t), y1(t), and z1(t) have the same number of elements
plot3(x1,y1,z1,'k','linewidth',1) % is used for 3-D plot
grid on
xlabel('x'); ylabel('y'); zlabel('z')

t=0:0.1:6*pi;
x1=sqrt(t).*sin(2*t);
y1=sqrt(t).*cos(2*t);

1
 z1=0.5*t;
x2=sqrt(t).*sin(4*t);
y2=sqrt(t).*cos(4*t);
z2=t;
x3=sin(t);
y3=cos(t);
z3=2*t;
plot3(x1,y1,z1,'k',x2,y2,z2,'r',x3,y3,z3,'m')
grid on
xlabel('x'); ylabel('y'); zlabel('z')

Mesh and Surface Plots

Mesh and surface plots visualize functions in 3-D. Here x and y independent variables.

• Mesh Plot

x=-8:0.5:8;
y=-10:0.5:10;
% Create a grid in the xy-plane of the cartesian coordinates:
[X,Y]=meshgrid(x,y) % The matrix X has the x coordinates of the grid points;
the matrix Y has the y coordinates of the grid points.
Y = 41×33
-10.0000 -10.0000 -10.0000 -10.0000 -10.0000 -10.0000 -10.0000 -10.0000

2
 -9.5000 -9.5000 -9.5000 -9.5000 -9.5000 -9.5000 -9.5000 -9.5000
-9.0000 -9.0000 -9.0000 -9.0000 -9.0000 -9.0000 -9.0000 -9.0000
-8.5000 -8.5000 -8.5000 -8.5000 -8.5000 -8.5000 -8.5000 -8.5000
-8.0000 -8.0000 -8.0000 -8.0000 -8.0000 -8.0000 -8.0000 -8.0000
-7.5000 -7.5000 -7.5000 -7.5000 -7.5000 -7.5000 -7.5000 -7.5000
-7.0000 -7.0000 -7.0000 -7.0000 -7.0000 -7.0000 -7.0000 -7.0000
-6.5000 -6.5000 -6.5000 -6.5000 -6.5000 -6.5000 -6.5000 -6.5000
-6.0000 -6.0000 -6.0000 -6.0000 -6.0000 -6.0000 -6.0000 -6.0000
-5.5000 -5.5000 -5.5000 -5.5000 -5.5000 -5.5000 -5.5000 -5.5000

Consider: over the domain and (as defined above).

Z=sqrt(X.^2 + Y.^2);
mesh(X,Y,Z)
colorbar
xlabel('x'); ylabel('y'); zlabel('z')

• Surface Plot

x=-8:0.5:8;
y=-10:0.5:10;
[X,Y]=meshgrid(x,y);
Z=sqrt(X.^2 + Y.^2);
surf(X,Y,Z)
colorbar

3
xlabel('x'); ylabel('y'); zlabel('z')

• Mesh Curtain

x=-8:0.5:8;
y=-10:0.5:10;
[X,Y]=meshgrid(x,y);
Z=sqrt(X.^2 + Y.^2);
meshz(X,Y,Z) % It surrounds the mesh with a curtain.
colorbar
xlabel('x'); ylabel('y'); zlabel('z')

4
 • Mesh and Contour

x=-8:0.5:8;
y=-10:0.5:10;
[X,Y]=meshgrid(x,y);
Z=sqrt(X.^2 + Y.^2);
meshc(X,Y,Z) % Below the mesh, a contour plot is created.
colorbar
xlabel('x'); ylabel('y'); zlabel('z')

5
 • Surface and Contour

x=-8:0.5:8;
y=-10:0.5:10;
[X,Y]=meshgrid(x,y);
Z=sqrt(X.^2 + Y.^2);
surfc(X,Y,Z) % Below the surface, a contour plot is created.
colorbar
xlabel('x'); ylabel('y'); zlabel('z')

6
7
 • Waterfall

x=-8:0.5:8;
y=-10:0.5:10;
[X,Y]=meshgrid(x,y);
Z=sqrt(X.^2 + Y.^2);
waterfall(X,Y,Z) % generates a mesh only in one direction.
colorbar
xlabel('x'); ylabel('y'); zlabel('z')

8
 • 3-D Contour

x=-8:0.5:8;
y=-10:0.5:10;
[X,Y]=meshgrid(x,y);
Z=sqrt(X.^2 + Y.^2);
contour3(X,Y,Z,10) % displays ''3-D surfaces by plotting z-slides on a 2-D
surface''. 10 is the number of contour levels, which is optional.
colorbar
xlabel('x'); ylabel('y'); zlabel('z')

9
 • 2-D Contour

x=-8:0.5:8;
y=-10:0.5:10;
[X,Y]=meshgrid(x,y);
Z=sqrt(X.^2 + Y.^2);
contour(X,Y,Z,10) % gives the projection of the contour levels on the xy-
plane(In this case, it is 10.)
colorbar
xlabel('x'); ylabel('y'); zlabel('z')

10
Plots with Special Graphics
• Plotting Sphere

% Make unit sphere

[x,y,z] = sphere; % Here x,y,z are the coordinates of the unit sphere.
surf(x,y,z) % draws the unit sphere

11
% Scale to desire radius.
r = 2;
% We introduce new coordinates X,Y,Z of a sphere with a radius of 2 by
multiplying the coordinates of the unit sphere.
X = x * r;
Y = y * r;
Z = z * r;
surf(X,Y,Z)
axis square % We introduce an axis style

12
 • Plotting Cylinder

% Make a cylinder with a radius of 1

[x,y,z]=cylinder; % Here x,y,z are the coordinates of the cylinder with a
radius of 1.
surf(x,y,z) % draws the cylinder with a radiues of 1

13
r = 4;
[x,y,z] = cylinder(r); % x,y,z are the new coordinates of the cylinder with
a radius of r=4.
h = 20;
Z = z*h; % We change the height of the cylinder.
surf(x,y,Z)

14
t = 0:pi/10:2*pi;
r = 2 + cos(t);
[x,y,z] = cylinder(r);
surf(x,y,z)
colorbar

15
 • 3-D Bar Plot

Z=[7 4 2 1];
bar3(Z) % draws 3-D bars of the elements of the vector Z.

16
Z= [1 2 5; 3 5 8; 5 8 9; 4 7 10];
bar3(Z) % draws 3-D bars of the elements of the matrix Z.

17
y = [1950 1960 1970 1980 1990]; % on the y-axis
z = [16 8 4 2 1]; % on the z-axis
bar3(y,z) % draws 3-D bars of the elements of the vector z at the values of
y.

18
 • 3-D Stem Plot

t=0:0.2:10;
x=t;
y=cos(t);
z=t*2;
stem3(x,y,z,'fill')
grid on
xlabel('x')
ylabel ('y')
zlabel('z')

19
 • 3-D Scatter Plot

t=0:0.2:10;
x=t;
y=cos(t);
z=t*2;
scatter3(x,y,z,'fill')
grid on
xlabel('x')
ylabel ('y')
zlabel('z')

20
 • 3-D Pie Plot

X=[20 18 10 5 32];
explode=[0 0 0 0 1]; % 1 separates the corresponding slice in X (in this
case the last slice) from the center.
pie3(X,explode)

21
Polar coordinates grid in the xy-plane:
t=(0:5:360)*pi/180; % angle of theta
r=0:0.1:2; % radius of r
% Create a grid of values of theta and radius of the polar coordinates:
[T,R]=meshgrid(t,r)
T = 21×73
0 0.0873 0.1745 0.2618 0.3491 0.4363 0.5236 0.6109
0 0.0873 0.1745 0.2618 0.3491 0.4363 0.5236 0.6109
0 0.0873 0.1745 0.2618 0.3491 0.4363 0.5236 0.6109
0 0.0873 0.1745 0.2618 0.3491 0.4363 0.5236 0.6109
0 0.0873 0.1745 0.2618 0.3491 0.4363 0.5236 0.6109
0 0.0873 0.1745 0.2618 0.3491 0.4363 0.5236 0.6109
0 0.0873 0.1745 0.2618 0.3491 0.4363 0.5236 0.6109
0 0.0873 0.1745 0.2618 0.3491 0.4363 0.5236 0.6109
0 0.0873 0.1745 0.2618 0.3491 0.4363 0.5236 0.6109
0 0.0873 0.1745 0.2618 0.3491 0.4363 0.5236 0.6109

R = 21×73
0 0 0 0 0 0 0 0
0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000 0.1000
0.2000 0.2000 0.2000 0.2000 0.2000 0.2000 0.2000 0.2000
0.3000 0.3000 0.3000 0.3000 0.3000 0.3000 0.3000 0.3000
0.4000 0.4000 0.4000 0.4000 0.4000 0.4000 0.4000 0.4000
0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000 0.5000
0.6000 0.6000 0.6000 0.6000 0.6000 0.6000 0.6000 0.6000
0.7000 0.7000 0.7000 0.7000 0.7000 0.7000 0.7000 0.7000
0.8000 0.8000 0.8000 0.8000 0.8000 0.8000 0.8000 0.8000
0.9000 0.9000 0.9000 0.9000 0.9000 0.9000 0.9000 0.9000

22
 Z=R.*T; % at each point of the grid
[X,Y]=pol2cart(T,R); % converts the polar coordinates grid to the Cartesian
coordinates grid
mesh(X,Y,Z)

On the view Command

23
x=-8:0.5:8;
y=-10:0.5:10;
[X,Y]=meshgrid(x,y);
Z=sqrt(X.^2 + Y.^2);
surf(X,Y,Z)
view(60,30) % changes the view of the plot
colorbar
xlabel('x'); ylabel('y'); zlabel('z')

24
25
SYBMOLIC MATH
A symbolic variable is a quantity that is represented by a symbol, such as a letter, rather than a numerical value.

Creating Symbolic Variables

x=3/7 % x stores 3/7 numerically
x = 0.4286
xs=sym(3/7) % xs stores 3/7 symbolically
xs =

syms x % is used to generate a symbolic variable 'x' without assigning it a

value
syms a b c y % is used to generate symbolic variables 'a b c y' without
assigning them values
x^2+3

ans =

a*sqrt(x)+c

ans =

2*b*exp(y)+x*c

ans =

Creating Symbolic Expressions

z=a*x^2+b*x+c % z is a symbolic expression

z =

g=3*x^2+2*x-1 % g is a symbolic expression

g =

f=b*exp(y)+x*c % f is a symbolic expression

f =

h=z^2+3*x*f+1 % h is a symbolic expression where the other symbolic

expressions, f and z, are used.

h =

• Substituting Values

1
v = subs(expression,variable,value)

subs(g,x,2) % substitute 2 for x in the expression g

ans =

subs(g,x,[2 3]) % substitute 2 and 3 for x in the expression g, separately

ans =

subs(z,b,3) % substitute 3 for b in the expression z

ans =

subs(f,[b c],[1 4]) % substitute 1 for b; substitute 4 for c in the

expression f

ans =

Manipulation of Symbolic Expressions

• On the command 'collect'

X=(2*x^2+x-3)*(3*x-exp(x))

X =

collect(X) % is used for ordering the terms in decreasing order of power of

x.

ans =

XY=(x^2-3*y^2-y*x+3*x)*(x*y^2-1)

XY =

collect(XY) % is used for ordering the terms in decreasing order of power of

x.

ans =

collect(XY,y) % is used for ordering the terms in decreasing order of power

of y.

ans =

• On the command 'expand'

expand(X)

ans =

2
 expand(cos(x+y))

ans =

• On the command 'factor'

F=x^3+4*x^2-11*x-30

F =

factor(F) % converts the polynomial expression F into a product of lower

degree polynomials.

ans =

Answer:

factor(2*x^2+x-3) % converts the polynomial expression into a product of

lower degree polynomials.

ans =

Answer:

• On the command 'simplify'

G=x*(x+2)^2-4*x

G =

simplify(G) % gives the simpler form of the expression G

ans =

simplify(cos(x)*cos(y)-sin(x)*sin(y)) % gives the simpler form of the

trigonometric expression

ans =

ans =

simplify((x^2+5*x+6)/(x+2)) % gives the simpler form of the expression

SOLVING ALGEBRAIC EQUATIONS

• Solving a Single Equation

syms a b c x y z
eq=(x/7)-3==5 % we define a symbolic equation 'eq'
eq =

3
xsol=solve(eq,x) % solves the symbolic equation 'eq' for x

xsol =

eq2=x^2-x-6==0 % we define a symbolic equation 'eq2'

eq2 =

xsol2=solve(eq2,x) % solves the symbolic equation 'eq2' for x

xsol2 =

T=a*x^2+5*b*x+20==0 % we define a symbolic equation 'T'

T =

Tx=solve(T,x) % solves the symbolic equation 'T' for x

Tx =

Ta=solve(T,a) % solves the symbolic equation 'T' for a

Ta =

• Solving a System of Equations

S=5*x+6*y+8*a==0 % we define a symbolic equation 'S'

S =

P=3*y-2*x==4*a % we define a symbolic equation 'P'

P =

[xSol,ySol]=solve([S P],[x y]) % solves the system of equations for x and y

xSol =

4
 ySol =

[aSol,ySol]=solve([S P],[a y]) % solves the system of equations for a and y

aSol =

ySol =

For different uses of 'solve' command, study the section 11.3 in your textbook !

Calculus with Symbolic Math

• Differentiation

syms x
y=(3+x)^3-exp(3*x)

y =

dy=diff(y,x) % gives the 1st derivative of y with respect to x

dy =

dy2=diff(y,x,2) % gives the 2nd derivative of y with respect to x

dy2 =

syms y
z=y*cos(3*x)-10*y^2*x

z =

dz=diff(z,y) % gives the 1st derivative of z with respect to y

dz =

dz=diff(z,x) % gives the 1st derivative of z with respect to x

dz =

syms n

5
f=x^n

f =

df=diff(f,x)

df =

% syms k
% dg=diff(y^k,y) % Warning: 'Unrecognized function or variable 'k''
% syms l
% g=3*l-a*x-3*y % Warning: 'Unrecognized function or variable 'l''

• Integration

y=(3+x)^3-exp(3*x)

y =

yint=int(y,x) % gives the integral of y with respect to x

yint =

yint03=int(y,x,[0 3]) % gives the integral of y with respect to x between

the end points 0 and 3
yint03 =

syms t
z=t*cos(3*x)-10*t^2*x

z =

zint=int(z,x) % gives the integral of z with respect to x

zint =

zint=int(z,t) % gives the integral of z with respect to t

zint =

zint=int(z,x,[1 4]) % gives the integral of z with respect to x between the

end points 1 and 4

6
 zint =

int(2*x,x)

ans =

Solving an Ordinary Differential Equation

• On first-order ODE

Consider:

syms y(t) % We define a symbolic function 'y(t)'

deqn=diff(y,t)+4*y==60 % We define a symbolic differential equation 'deqn'

ysol=dsolve(deqn) % gives general solution to 'deqn'

ysol =

can be obtained if initial condition is given: Consider

ic=y(0)==5 % initial condition

ic =

yp=dsolve(deqn,ic) % gives particular solution to 'deqn'

yp =

Another Method
dsolve('Dy+4*y=60') % Here Dy means dy/dt
Warning: Support for character vector or string inputs will be removed in a future release.
Instead, use syms to declare variables and replace inputs such as dsolve('Dy = -3*y') with syms
y(t); dsolve(diff(y,t) == -3*y).

ans =

dsolve('Dy+4*y=60','y(0)=5')
Warning: Support for character vector or string inputs will be removed in a future release.
Instead, use syms to declare variables and replace inputs such as dsolve('Dy = -3*y') with syms
y(t); dsolve(diff(y,t) == -3*y).

7
 ans =

• On higher-order ODE

Consider:

deqn2=diff(y,t,2)-2*diff(y,t)+2*y==0 % We define symbolic differential

equation 'deqn2'.
deqn2(t) =

ysol2=dsolve(deqn2) % gives general solution to 'deqn2'

ysol2 =

and can be obtained if initial conditions are given: Consider

ic1=y(0)==1 % initial condition 1

ic1 =

Ay=diff(y)
Ay(t) =

ic2=Ay(0)==0 % initial condition 2

ic2 =

yp2=dsolve(deqn2,ic1,ic2) % gives particular solution to 'deqn2'

yp2 =

Another Method
dsolve('D2y-2*Dy+2*y=0') % Here D2y=d^2y/dt^2

8
 Warning: Support for character vector or string inputs will be removed in a future release.
Instead, use syms to declare variables and replace inputs such as dsolve('Dy = -3*y') with syms
y(t); dsolve(diff(y,t) == -3*y).

ans =

dsolve('D2y-2*Dy+2*y=0','y(0)=1','Dy(0)=0')
Warning: Support for character vector or string inputs will be removed in a future release.
Instead, use syms to declare variables and replace inputs such as dsolve('Dy = -3*y') with syms
y(t); dsolve(diff(y,t) == -3*y).

ans =

Plotting Symbolic Functions and Equations

syms x y

h=x^2+3*x; % 'h' is a symbolic function

fplot(h) % is used to plot the symbolic function 'h'
xlabel('x')
ylabel('h')

9
fimplicit(h==10,[-10,10]) % is used to plot the equation for 'h=10'

u=x^2+y^2;
fimplicit(u==25) % is used to plot the equation for 'u=25'
hold on
fimplicit(u==2) % is used to plot the equation for 'u=2'
hold off

10
fcontour(u) % is used to plot the equations for many values of u

11
fcontour(u,"LevelList",[1 4 9 16]) % is used to plot the equations for some
specified values of u

12
Another Method
g=3*x^2+5;
ezplot(g) % is used to plot g(x) for one variable x

13
g=3*x^2+5;
ezplot(g,[-3,3]) % The domain of the variable x is btw -3 and 3

14
g=3*x^2+5;
ezplot(g,[-3,3,10,30]) % The domain of the variable x is btw -3 and 3; the
domain of the variable g is btw 10 and 30

15
f=x^2+y^2-5;
ezplot(f) % is used to plot f(x,y)=0 for two variables x and y.

16
f=x^2+y^2-5;
ezplot(f,[-4,4,-4,4])

17
syms t
x=sin(2*t)-5;
y=cos(4*t)+3;
ezplot(x,y)

18
19