INTRODUCTION TO MATLAB

1. I NSTALLATION

The National University of Singapore has a Total Academic Headcount licence for MATLAB.
Students may use it for academic, research, and learning. The license allows students to install
MATLAB on personally-owned computers.
1. If you are NOT using NUS network, you are required to use nVPN:
◦ https://webvpn.nus.edu.sg/
Please sign in with your NUSNET ID in the format: nusstu\nusnetid and password.
If you are using NUS network, please proceed to Step 2.

2. Click the link https://sm05.stf.nus.edu.sg/studentmatlab/ to download MATLAB. It

is required to sign in with your NUSNET ID in the format: nusstu\nusnetid and password.
Then follow the instructions.

3. Create a MathWorks account using your NUS email address.

(i) https://www.mathworks.com/mwaccount/register
¦ Email address: Your NUS email account (e.g., e0012345@u.nus.edu)
¦ Location: Singapore
¦ How will you use MathWorks software? Student use
¦ Are you at least 13 years or older? Yes
(ii) You will receive an email from service@mathworks.com with title “Verify Email Ad-
dress”. Click the link in the email to verify your account.
(iii) Finish creating your profile. Then you should be able to see the following information:
◦ Your account has been created and license 40707750 has been associated with your
account.

4. Click the Download bottom to download and run the installer.

(i) When prompted, log in with your MathWorks Account (your NUS email account).
(ii) Select your licence (40707750, Student, Academic — Total Headcount).
(iii) Choose installation folder.
(iv) Select products to install.
1
 INTRODUCTION TO MATLAB 2

2. B ASIC O PERATIONS

First of all, we learn some basic operations in MATLAB.

MATLAB environment behaves like a super-complex calculator. You can enter the com-
mands at the >> command prompt. The answer appears by pressing Enter.
We can use the following arithmetic operators: addition + , subtraction - , multiplication
* , division / , exponentiation  .
For example, to add two numbers 123 and 321, we simply type
>> 123 + 321
and then press Enter:
ans = 444
Similarly, we can apply other operators to use MATLAB as an ordinary calculator:
>> 123 - 321
ans = -198
>> 123 * 321
ans = 39483
>> 123 / 321
ans = 0.3832
>> 123  3
ans = 1860867
You may add a semicolon ; at the end of the statement; then MATLAB will hide the output.
For example,
>> a = 3;
>> a  2
ans = 9
Note that the symbol a is now defined as 3. We may remove it from the memory by using
clear a or remove all variables from the memory by clear . If we want to clear the command
window, we can use clc .

3. B ASIC F UNCTIONS

The function sqrt(x) computes the principle square root of the number x. For example,
>> sqrt(2)
 INTRODUCTION TO MATLAB 3

ans = 1.4142
By default, MATLAB displays four decimal digits to its answers. But we can change the for-
mat for numeric display. (The percent symbol % is used for indicating a comment line.)
>> format long % 16 decimal digits
>> sqrt(2)
ans = 1.414213562373095
>> format rat % rational approximation
>> sqrt(2)
ans = 1393/985
>> format short % four decimal digits (default)
>> sqrt(2)
ans = 1.4142
We can also use the following functions in MATLAB:
(i) Trigonometric functions: sin(x) , cos(x) , tan(x) , cot(x) , sec(x) , csc(x) .
(ii) Inverse trigonometric functions: asin(x) , acos(x) , atan(x) , acot(x) .
(iii) Exponential and logarithmic functions: exp(x) , log(x) (base e), log10(x) (base 10).

For trigonometric and inverse trigonometric functions, the angles are measured in radian (π
radian = 180◦ ), and the constant π = 3.1415926 · · · is defined by pi ( MATLAB is case-sensitive).
>> sin(pi/4)
ans = 0.7071
>> atan(1)
ans = 0.7854

4. S OLVING B ASIC A LGEBRAIC E QUATIONS

The command solve is used for solving algebraic equations. We shall

(i) Use syms to declare the variables.

(ii) Use solve(equations, variables) to solve the specific equations.

Note that = is used to assign a value to a variable; and == shall be used for equality.

4.1. Single Variable Equations. For example, x 2 + x − 1 = 0.

>> syms x;
>> solve(x2 + x - 1 == 0, x)
 INTRODUCTION TO MATLAB 4

ans = -5(1/2)/2 - 1/2

ans = 5(1/2)/2 - 1/2
We may use vpa(x) or vpa(x,d) to evaluate each element of the symbolic input to d
digits (the default is 32 digits).
>> vpa(ans,10)
ans = -1.618033989
ans = 0.6180339887
Sometimes the exact solutions cannot be specifically displayed. For example,

x 4 − 7x 3 + 3x 2 − 5x + 9 = 0.

>> solve(x4 - 7*x3 + 3*x2 - 5*x + 9 == 0, x)

ans = root(z4 - 7*z3 + 3*z2 - 5*z + 9, z, 1)
ans = root(z4 - 7*z3 + 3*z2 - 5*z + 9, z, 2)
ans = root(z4 - 7*z3 + 3*z2 - 5*z + 9, z, 3)
ans = root(z4 - 7*z3 + 3*z2 - 5*z + 9, z, 4)
Then evaluate using floating points to obtain the decimal expressions of the four roots:
>> vpa(ans,10)
ans = 1.059780463
ans = - 0.3450883978 - 1.077836295i
ans = - 0.3450883978 + 1.077836295i
ans = 6.630396332
Here the symbol i refers to the imaginary unit i where i = e πi /2 such that i 2 = −1.
Alternatively, one may use vpasolve to get the decimal expression directly:
>> vpasolve(x4 - 7*x3 + 3*x2 - 5*x + 9 == 0, x)
ans = 1.0597804633025896291682772499885
ans = 6.630396332390718431485053218985
ans = - 0.34508839784665403032666523448675 - 1.0778362954630176596831109269793i
ans = - 0.34508839784665403032666523448675 + 1.0778362954630176596831109269793i
MATLAB can also be used to solve symbolic equations. For example, ax 2 + bx + c = 0.
>> syms a b c x;
>> solve(a*x2 + b*x + c == 0, x)
ans = -(b + (b2 - 4*a*c)(1/2))/(2*a)
 INTRODUCTION TO MATLAB 5

ans = -(b - (b2 - 4*a*c)(1/2))/(2*a)

4.2. Multi-Variable Equations. MATLAB can also solve multi-variable equations, for which
the equations shall be put together in square brackets [ ] , so do the variables. Moreover, since
the output consists multiple variables, we shall assign the solution as a vector. For example,

3x + y = 10, x + y = 20.

>> syms x y;
>> [Sx,Sy] = solve([3*x+y==10, x+y==20], [x,y])
Sx = -5
Sy = 25
The solution of the a general linear system in variables x and y
(
ax + b y = e
cx + d y = f
can be obtained as follows:
>> syms a b c d e f x y;
>> [Sx,Sy] = [Sx, Sy] = solve([a*x+b*y == e, c*x+d*y == f], [x,y])
Sx = -(b*f - d*e)/(a*d - b*c)
Sy = (a*f - c*e)/(a*d - b*c)

5. S UMS

If we need to find the sum of a sequence, use the command symsum . The format is
>> symsum(expression, variable, [min, max])
For example, for the arithmetic sequence {a n } given by a n = a + (n − 1)d , the sum of its first
n terms is (2a + (n − 1)d )n/2.
>> syms a d k n;
>> symsum(a+(k-1)*d, k, [1, n])
ans = a*n - d*(n - (n*(n + 1))/2)
>> simplify(ans) % simplify the previous answer
ans = (n*(2*a - d + d*n))/2
For the geometric sequence {a n } given by a n = ar n−1 , the sum of its first n terms is
(
a(r n − 1)/(r − 1) if r 6= 1,
Sn =
an if r = 1.
 INTRODUCTION TO MATLAB 6

>> syms a r k n;
>> symsum(a*r(k-1), k, [1, n])
ans = piecewise(r == 1, a*n, r ∼= 1, (a*(rn - 1))/(r - 1))
Moreover, we can find the sum to infinity, which is defined by inf in MATLAB. Recall that
n−1
for the geometric sequence {a n } given by a n = ar , the sum to infinity is
(
a/(1 − r ) if |r | < 1,
S∞ =
does not exist if |r | ≥ 1.

>> syms a r k;
>> assume(abs(r)<1);
>> symsum(a*r(k-1), k, [1, inf])
ans = -a/(r - 1)

6. V ECTOR

A vector a i + b j + c k or (a, b, c) can be defined as [a, b, c] or simply [a b c] .

The addition, subtraction and multiplication with numbers can be evaluated using + , -
and + respectively.
For example, let u = (1, 2, 3) and v = (4, 5, 6).
>> u = [1 2 3];
>> v = [4 5 6];
>> u + v
ans = 5 7 9
>> u - v
ans = -3 -3 -3
>> 3*u
ans = 3 6 9
The dot product u • v of two vectors u and v is defined by dot(u,v) .
>> dot(u,v)
ans = 32
p
We can use the formula |v | = v • v to find the norm of v :
>> sqrt(dot(v,v))
ans = 8.7750
 INTRODUCTION TO MATLAB 7

Alternatively, MATLAB provides a command norm for the norm of a vector.

>> norm(v)
ans = 8.7750
The cross product u × v is defined by cross(u,v) .
>> cross(u,v)
ans = -3 6 -3

7. F UNCTION

7.1. Standard Function. To define a function, we shall (i) give the name of the function, (ii) use
@ to declare the name of the variable, (iii) provide the expression of the function. For example,
define f (x) = x 2 :
>> f = @(x) x2;
Then f (2) can be evaluated by
>> f(2);
ans = 4
Another example: g (x) = sin(x 3 )/x 2
>> g = @(x) sin(x3)/x2;

Multi-variable functions can be defined similarly by declaring more variables. For example,
p
h(x, y) = x 2 + y 2 :
>> h = @(x,y) sqrt(x2 + y2);
To evaluate h(5, 12):
>> h(5,12)
ans = 13

7.2. Piecewise Function. The absolute value function f (x) = |x| is a piecewise function:
(
x if x ≥ 0,
f (x) =
−x if x < 0,

which can be defined in MATLAB:

>> syms x
>> f = piecewise(x>=0, x, x<0, -x);
The general format is to use the piecewise as
>> f = piecewise(condition 1, value 1, ..., condition N, value N);
 INTRODUCTION TO MATLAB 8

In the example above, the first condition is x>=0 and this means “x is greater than or equal
to 0”. The value x is the expression of the function when the condition x>=0 is satisfied. It
defines f (x) = x when x ≥ 0. The second condition x<0 and the second value -x defines
f (x) = −x when x < 0.
In order to evaluate the value of the function at given point, for example, f (−2), we use
>> subs(f, x, -2)
ans = 2

8. C URVE P LOTTING

8.1. Standard Function. Once a function is defined, fplot can be used to plot its graph over
a specified interval in the form [xmin, xmin] .
For example, we plot f (x) = x 2 on the interval (−100, 100):
>> f = @(x) x2;
>> fplot(f, [-100,100]);

The command can be used even if the function is undefined somewhere on the specific in-
terval. For example, g (x) = 1/x on (−1, 1) is defined at x = 0.
>> g = @(x) 1/x;
>> fplot(g, [-1,1]);
 INTRODUCTION TO MATLAB 9

8.2. Multiple Curves. In order to plot more curves in the same coordinate system, we use the
commands on and off . For example,

f (x) = sin x and g (x) = cos t , −2π < x < 2π.

(We can use Shift + Enter to start a new line without breaking the command.)
>> f = @(x) sin(x);
>> g = @(x) cos(x);
>> fplot(f, [-2*pi, 2*pi]);
>> hold on
>> fplot(g, [-2*pi, 2*pi]);
>> hold off

We can draw multiple graphs on the same plot. MATLAB provides some basic colour options:
white w , black b , blue b , red r , cyan c , green g , magenta m , yellow y to distinguish
 INTRODUCTION TO MATLAB 10

different graphs. For example,

f (x) = x 3 − x + 1,

g (x) = x 4 − 3x 2 + x,

h(x) = x 5 + 0.3x 4 − 2.8x 3 − 0.3x 2 + 1.8x.

>> f = @(x) x3 - x + 1;

>> g = @(x) x4 - 3*x2 + x;
>> h = @(x) x5 + 0.3*x4 - 2.8*x3 - 0.3*x2 + 1.8*x;
>> fplot(f, [-2,2], 'r');
>> hold on
>> fplot(g, [-2,2], 'b');
>> fplot(h, [-2,2], 'g');
>> hold off;

8.3. Parametric Equations. A curve may be defined using a pair of parametric equations:

x = x(t ) and y = y(t ).

In order to use fplot , we shall first define the functions for the x- and y-coordinates.
Recall that a unit circle can be parameterized by

x = cos t , y = sin t , 0 ≤ t ≤ 2π.

>> f = @(t) cos(t);

>> g = @(t) sin(t);
>> fplot(f,g, [0,2*pi]);
 INTRODUCTION TO MATLAB 11

MATLAB can draw a butterfly:

>> f = @(t) sin(t)*(exp(cos(t))-2*cos(4*t) - (sin(t/12))5);
>> g = @(t) cos(t)*(exp(cos(t))-2*cos(4*t) - (sin(t/12))5);
>> fplot(f,g, [0,12*pi]);

8.4. Piecewise Function. There are two ways to plot a piecewise function.
Once a piecewise function is already defined, fplot can directly plot its graph:

(
ex if x < 0,
f (x) =
cos x if x ≥ 0.

>> syms x;
>> f = piecewise(x<0, exp(x), x>=0, cos(x));
>> fplot(f, [-3,3]);
 INTRODUCTION TO MATLAB 12

Alternatively, we can plot different branches on different intervals using the same colour.
>> f1 = @(x) exp(x);
>> f2 = @(x) cos(x);
>> fplot(f1, [-3,0], 'b');
>> hold on
>> fplot(f2, [0,3], 'b');
>> hold off;

8.5. Implicit Function. The command fimplicit can be used to plot graphs defined by im-
plicit functions, i.e., f (x, y) = 0. The format is:
>> fimplicit(function in two variables, [xmin, xmax, ymin, ymax])
For example, to plot x 3 + y 3 = 3x y, we shall first define the function f (x, y) = x 3 + y 3 − 3x y:
>> f = @(x,y) x3 + y3 - 3*x*y;
>> fimplicit(f, [-1.5, 2, -1.5, 2]);
 INTRODUCTION TO MATLAB 13

9. L IMITS

The limit of a function describes the behavior of a function near a point or at infinity. The
limit lim f (x) can be easily evaluated with the command limit :
x→a
>> limit(f(x), x, a)
x
For example, lim p :
x→0 1 + 3x − 1
>> syms x;
>> limit(x/(sqrt(1+3*x)-1), x, 0)
ans = 2/3
The one-sided limits can be evaluated by specifying the direction left or right . For ex-
ample, let the floor function bxc denote the greatest integer smaller or equal to x. Then

lim bxc = 1, lim bxcx = 2, and lim bxc does not exist.
x→2− x→2+ x→2

>> syms x;
>> limit(floor(x), x, 2, 'left')
ans = 1
>> limit(floor(x), x, 2, 'right')
ans = 2
>> limit(floor(x), x, 2)
ans = limit(floor(x), x, 2)
³ a ´x
We can find the limit at infinity inf or negative infinity -inf . For example, lim 1 + :
x→∞ x
>> syms a x;
>> limit((1+a/x)x, x, inf)
 INTRODUCTION TO MATLAB 14

ans = exp(a)
Sometimes the limit depends on the value of the parameters. For example,

 0 if a < 0,


lim 2ax = 1 if a = 0,
x→∞ 

 ∞ if a > 0.

We may use assume to declare the value of the parameter.

>> syms a x;
>> assume(a<0)
>> limit(2(a*x), x, inf)
ans = 0
>> assume(a==0)
>> limit(2(a*x), x, inf)
ans = 1
>> assume(a>0)
>> limit(2(a*x), x, inf)
ans = Inf

10. D ERIVATIVE

The derivative of a function f at point a is defined as the limit

f (x + h) − f (x)
f 0 (x) = lim .
h→0 h

For example, if f (x) = x 2 , then f 0 (x) = 2x:

>> syms x h;
>> f = @(x) x2;
>> limit((f(x+h)-f(x))/h, h, 0)
ans = 2*x
On the other hand, MATLAB includes the command diff for differentiation. For example,
2
find the derivative of f (x) = cos(x ) with respect to x:
>> syms x;
>> f = @(x) cos(x2);
>> diff(f, x)
ans = -2*x*sin(x2)
 INTRODUCTION TO MATLAB 15
p
x − 3x x
In order to differentiate g (x) = p with respect to x:
x
>> syms x;
>> g = @(x) (x-3*x*sqrt(x))/sqrt(x);
>> diff(g, x)
ans = - (x - 3*x(3/2))/(2*x(3/2)) - ((9*x(1/2))/2 - 1)/x(1/2)
The answer seems complicated. Fortunately MATLAB provides the command simplify which
may help in this situation:
>> simplify(ans)
ans = -(6*x(1/2) - 1)/(2*x(1/2))
Suppose we are looking for f (4) (x), the 4th order derivative of f (x) with respect to 4, of course
we may use
>> diff(diff(diff(diff(f, x), x), x), x)
ans = 16*x4*cos(x2) - 12*cos(x2) + 48*x2*sin(x2)
However, it may be a bit cumbersome. The following shorter commands have the same effect:
>> diff(f, x,x,x,x)
ans = 16*x4*cos(x2) - 12*cos(x2) + 48*x2*sin(x2)
>> diff(f, x, 4)
ans = 16*x4*cos(x2) - 12*cos(x2) + 48*x2*sin(x2)
For example, we can find the local extreme values of f (x) = x 3 − 3x 2 + x − 2.
>> syms x;
>> f = @(x) x3 - 3*x2 + x -2;
>> g = diff(f, x); % define g(x) to be f'(x)
>> h = diff(f, x, x); % define h(x) to be f(x)
(i) Solve f 0 (x) = 0 to obtain the stationary points.
>> solve(g(x) == 0, x)
ans = 1 - 6(1/2)/3
ans = 6(1/2)/3 + 1
(ii) Check the sign of f 00 (x) at the stationary points. (Note that you can use Copy and Paste
in the scripts.)
>> subs(h, x, 1-6(1/2)/3)
ans = -22063042185692343/4503599627370496
 INTRODUCTION TO MATLAB 16

>> subs(h, x, 6(1/2)/3+1)

ans = 5515760546423085/1125899906842624
p p
Hence, f has a local maximum at 1 − 6/3 and a local minimum at 1 + 6/3.
We can plot f (x), f 0 (x) and f 00 (x) in the same coordinate system.
>> fplot(f, [-2,5], 'r');
>> hold on
>> fplot(g, [-2,5], 'b');
>> fplot(h, [-2,5], 'g');
>> hold off

11. I NTEGRAL

11.1. Indefinite Integral. For a function f (x), its indefinite integral is a function F (x) such that
F 0 (x) = f (x), denoted by
Z
F (x) = f (x) d x.

All indefinite integrals differ by a constant only. So we also use

Z
F (x) = f (x) d x +C

to represent the entire family of indefinite integrals of f (x), where C is an arbitrary constant.
In MATLAB, we can use the command int to find an indefinite integral.
1 1 x
Z
For example, 2 2
d x = tan−1 x + 2
+C .
(1 + x ) 2 2(x + 1)
>> syms x;
>> f = @(x) 1/(1+x2)2;
>> int(f, x)
ans = atan(x)/2 + x/(2*(x2 + 1))
 INTRODUCTION TO MATLAB 17

Note that the constant C is dropped in the answer. Another example:

5 15 3 1
Z
sin6 x d x = x − sin 2x + sin 4x − sin 6x +C .
16 64 64 192
>> syms x;
>> g = @(x) sin(x)6;
>> int(g, x)
ans = (5*x)/16 - (15*sin(2*x))/64 + (3*sin(4*x))/64 - sin(6*x)/192
Z b
11.2. Definite Integral. The definite integral f (x) d x represents the net area bounded be-
a
tween the graph of y = f (x) and the x-axis from x = a to x = b. The command int can also be
used to find definite integral by specifying the lower limit a and upper limit b in the form
>> int(expression, variable, [lower limit, upper limit])
1 1 π 1
Z
For example, 2 2
dx = + .
0 (1 + x ) 8 4
>> syms x;
>> f = @(x) 1/(1+x2)2;
>> int(f, x, [0, 1])
ans = pi/8 + 1/4
Another example: Z π 5
sin6 x d x = π.
0 16
>> syms x;
>> g = @(x) sin(x)6;
>> int(g, x, [0, pi])
ans = (5*pi)/16
We can verify the fundamental theorem of calculus:
Z x
d
f (t ) d x = f (x).
dx a
>> syms a t x f(t); % use f(t) to declare that f is a function in t
>> diff(int(f(t), t, [a, x]), x)
ans = f(x)
It is also possible to integrate over an infinite interval by specifying the lower limit and/or the
upper limit as -inf or inf . For example,
Z ∞
2 p
e −x d x = π.
−∞
 INTRODUCTION TO MATLAB 18

>> int(exp(-x2), x, [-inf, inf])

ans = pi(1/2)

12. D IFFERENTIAL E QUATION

dy
12.1. General Solution. In order to solve an ordinary differential equation of the form =
dx
f (x, y), we shall first declare that y is a function in x, and then use dsolve .
dy
For example, = 1 + x + y. Note that the equal sign is represented by == .
dx
>> syms y(x); % declare that y is a function in x
>> dsolve(diff(y, x) == 1+x+y)
ans = C1*exp(x) - x -2
Here C1 refers to an undermined constant. So the general solution of the given equation is
y = C e x − x − 2, where C is an arbitrary constant.
Ordinary differential equations with higher order can be solved similarly.
d2y dy
For example, 2
+2 + 2y = e x .
dx dx
>> syms y(x);
>> dsolve(diff(y,x,x) + 2*diff(y,x) + 2*y == 0)
ans = C2*exp(-x)*cos(x) - C3*exp(-x)*sin(x)
The general solution is y = C 2 e −x sin x +C 3 e −x cos x, where C 2 ,C 3 are arbitrary constants.

12.2. Particular Solution. Suppose an initial condition is given. We can use

>> dsolve(ode, condition)
dy
For example, = 1 + x + y such that y = 1 at x = 0.
dx
>> syms y(x);
>> dsolve(diff(y, x) == 1+x+y, y(0)==1)
ans = 3*exp(x) - x-2
For higher order differential equations, we need more initial conditions. We use
>> dsolve(ode, [condition 1, ..., condition N])
The condition y 0 (a) = b shall be defined as Dy(a)=b if we set
>> Dy = diff(y,x);
d2y dy
For example, 2
+8 − 9y = 0, where y(1) = 2 and y 0 (1) = 0.
dx dx
>> syms y(x);
 INTRODUCTION TO MATLAB 19

>> Dy = diff(y,x);
>> dsolve(diff(y,x,x) + 8*diff(y,x) - 9*y == 0, [y(1) == 2, Dy(1) == 0])
ans = (exp(-9*x)*exp(9))/5 + (9*exp(-1)*exp(x))/5
>> simplify(ans)
ans = (9*exp(x - 1))/5 + exp(9 - 9*x)/5

12.3. System of Differential Equations. The command dsolve can also be used for a system
of differential equations. Note that the output consists multiple functions. So the output must
be in vector form:
>> [var 1, ..., var N] = dsolve([eqn1, ..., eqn M], [cond 1, ..., cond k])
For example, the following linear system consists of two functions x(t ) and y(t ) with two
equations and two initial conditions:
(
x 0 (t ) = 3x + 4y
where x(0) = 2 and y(0) = 3.
y 0 (t ) = −4x + 3y

>> syms x(t) y(t);

>> [Sx, Sy] = dsolve([diff(x,t) == 3*x + 4*y, diff(y,t) == -4*x + 3*y], [x(0)
== 2, y(0) == 3])
Sx = 2*cos(4*t)*exp(3*t) + 3*sin(4*t)*exp(3*t)
Sy = 3*cos(4*t)*exp(3*t) - 2*sin(4*t)*exp(3*t)

13. P LOTTING IN THE S PACE

13.1. Curve in the Space. A curve in the x y z-space can be parametrized as

r (t ) = x(t )i + y(t )j + z(t )k, a ≤ t ≤ b.

Plotting a space curve in the x y z-space is similar to plotting a parametrized curve in the x y-
plane. The only difference is to replace fplot by fplot3 . For example,

r (t ) = (e −t /10 sin 5t ) i + (e −t /10 cos 5t ) j + t k.

>> x = @(t) exp(-t/10)*sin(5*t);

>> y = @(t) exp(-t/10)*cos(5*t);
>> z = @(t) t;
>> fplot3(x,y,z, [-10,10])
 INTRODUCTION TO MATLAB 20

13.2. Level Curves. A surface is defined as a two-variable function z = f (x, y). Its level curve at
z = c is the intersection of the surface with the horizontal plane at z = c; that is, f (x, y) = c.
The command fcontour plots the level curve of a two-variable function.
>> fcontour(function, [xmin, xmax, ymin, ymax])
For example, f (x, y) = sin x + cos y:
>> f = @(x,y) sin(x) + cos(y);
>> fcontour(f, [-10,10, -10,10])

It is also possible plot the level curves of two surfaces using hold on and hold off .

13.3. Surface. A surface z = f (x, y) in two variables can be plotted in x y z-space using fsurf .
For example, f (x, y) = sin x + cos y.
>> f = @(x,y) sin(x) + cos(y);
>> fsurf(f, [-10,10, -10,10])
 INTRODUCTION TO MATLAB 21

13.4. Parametrized Surface. A surface can be parametrized by two parameters u and v:

r (u, v) = x(u, v) i + y(u, v) j + z(u, v) k.

We can still use fsurf to plot parametrized surface in the form

>> fsurf(x, y, z, [umin, umax, vmin, vmax])
For example,

r (u, v) = u cos v i + u sin v j + v k.

>> x = @(u,v) u*cos(v);

>> y = @(u,v) u*sin(v);
>> z = @(u,v) v;
>> fsurf(x,y,z, [-5,5, -5,5])
 INTRODUCTION TO MATLAB 22

13.5. Implicit Function. A surface can also be defined by an implicit function, i.e., f (x, y, z) = 0.
It is almost the same as plotting a curve defined by implicit function in the x y-plane. Simply
use fimplicit3 instead of fimplicit .
For example, plot the sphere x 2 + y 2 + z 2 = 4z and the paraboloid z = x 2 + y 2 .
>> f = @(x,y,z) x2 + y2 + z2 - 4*z;
>> g = @(x,y,z) x2 + y2 - z;
>> fimplicit3(f, [-3,3, -3,3, -1,5]);
>> hold on
>> fimplicit3(g, [-3,3, -3,3, -1,5]);
>> hold off

14. M ULTI -VARIABLE D ERIVATIVES

14.1. Partial Derivatives. For a multi-variable function, we can find its partial derivatives by
specifying the variables to be differentiated with respect to properly.
∂f ∂f ∂f
For example, for f (x, y, z) = x ln(x y 2 z 3 ), its partial derivatives , and are given by
∂x ∂y ∂z
>> syms x y z;
>> f = @(x,y,z) x*log(x*y2*z3);
>> diff(f, x)
ans = log(x*y2*z3)+1
>> diff(f, y)
ans = (2*x)/y
>> diff(f, z)
 INTRODUCTION TO MATLAB 23

ans = (3*x)/z
The command gradient provides the gradient vector ∇ f of a function f so that we can
obtain all the derivatives at once (the output is a column vector).
>> gradient(f, [x,y,z])
ans = log(x*y2*z3)+1
ans = (2*x)/y
ans = (3*x)/z
∂2 f ∂2 f
We may verify the mixed derivative theorem: = :
∂y ∂x ∂x ∂y
>> diff(f, x, y)
ans = 2/y
>> diff(f, y, x)
ans = 2/y
Recall that the directional derivative of f (x, y, z) along a vector v at (a, b, c) is given by

∇ f (a, b, c) • v /|v |.

For example, at (1, 2, 3), the directional derivative of f along v = (4, 5, 6) is

>> v = [4 5 6];
>> subs(gradient(f, [x,y,z]), {x,y,z}, {1,2,3});
>> dot(ans, v)/norm(v)
ans = (77(1/2)*(4*log(108) + 15))/77

14.2. Derivative of a Curve. For a curve C defined by

r (t ) = x(t )i + y(t )j + z(t )k,

its derivative vector is r 0 (t ) = x 0 (t )i + y 0 (t )j + z 0 (t )k.

One may define the x-, y- and z-components separately. Alternatively, we can define r as a
vector function in t . For example, r (t ) = (cos t )i + (sin t )j + t k.
>> r = @(t) [cos(t) sin(t) t];
Then its derivative vector is
>> syms t;
>> diff(r, t)
ans = [ -sin(t), cos(t), 1]
 INTRODUCTION TO MATLAB 24

15. I NTEGRAL FOR M ULTI -VARIABLE F UNCTIONS

15.1. Double Integral. Suppose that a region D in the x y-plane is parametrized by

a ≤ x ≤ b, α(x) ≤ y ≤ β(x).

Then for any function f (x, y), the double integral

Ï Z b Z β(x)
f (x, y) d A = f (x, y) d y d x,
D a α(x)

is converted to two integrals.

Ï Z 2Z 2x
2 2
2
Let D = {(x, y) | 0 ≤ x ≤ 2, x ≤ y ≤ 2x}, then (x + y ) d A = (x 2 + y 2 ) d y d x.
D 0 x2
>> syms x y;
>> f = @(x,y) x2 + y2;
Z 2x
We first find f (x, y) d y:
x2
>> int(f, y, [x2,2*x])
ans = -(x3*(x3 + 3*x - 14))/3
then integrate the result with respect to x on [0, 2]:
>> int(ans, x, [0,2])
ans = 216/35
p
Note that D can also be parametrized by D = {(x, y) | 0 ≤ y ≤ 4, y/2 ≤ x ≤ y}. So
Ï Z 4 Z py
2 2
(x + y ) d A = (x 2 + y 2 ) d x d y.
D 0 y/2

>> int(f, x, [y/2, sqrt(y)])

ans = y(3/2)*(y - (13*y(3/2))/24 + 1/3)
>> int(ans, y, [0,4])
ans = 216/35

15.2. Line Integral of a Function. Let C a curve parametrized by

r (t ) = x(t )i + y(t )j + z(t )k, a ≤ t ≤ b.

Its length is the line integral

Z Z b
1ds = |r 0 (t )| d t .
C a
For example, r = (cos t )i + (sin t )j + t k, 0 ≤ t ≤ 2π.
>> r = @(t) [cos(t) sin(t) t];
>> syms t;
 INTRODUCTION TO MATLAB 25

>> norm(diff(f,t))
ans = (abs(cos(t))2 + abs(sin(t))2 + 1)(1/2)
>> int(ans, [0, 2*pi])
ans = 2*pi*2(1/2)
In general, the line integral of a function f along C is
Z Z b
f ds = f (x(t ), y(t ), z(t )) |r 0 (t )| d t .
C a
Since there is no simple way to evaluate f (r (t )) directly, we may need to define all the compo-
nents of r as functions.
For example, r = (cos t )i + (sin t )j + t k, 0 ≤ t ≤ 2π, and f (x, y, z) = x 2 + y z.
>> x = @(t) cos(t);
>> y = @(t) sin(t);
>> z = @(t) t;
>> r = @(t) [x(t) y(t) z(t)]; % define the curve
>> f = @(x,y,z) x2 + y*z;
>> syms t;
>> dot(f(x(t),y(t),z(t)), norm(diff(r,t)));
>> int(ans, t, [0, 2*pi])
ans = -pi*2(1/2)

15.3. Line Integral of a Vector Field. A vector field F is a vector valued function in multi-
variables:
F (x, y, z) = P (x, y, z)i +Q(x, y, z)j + R(x, y, z)k.
Its line integral along the curve C defined by r (t ) = x(t )i + y(t )j + z(t )k, a ≤ t ≤ b, is
Z Z b
F • dr = F (x(t ), y(t ), z(t )) • r 0 (t ) d t .
C a

For example, r = (cos t )i + (sin t )j + t k, 0 ≤ t ≤ 2π, and F (x, y, z) = x 2 i + y 2 .

>> x = @(t) cos(t);
>> y = @(t) sin(t);
>> z = @(t) t;
>> r = @(t) [x(t) y(t) z(t)];
>> F = @(x,y,z) x2 + y2 + z2;
>> syms t;
 INTRODUCTION TO MATLAB 26

>> dot(F(x(t),y(t),z(t)), diff(r,t));

>> int(ans, t, [0, 2*pi])
ans = (8*pi3)/3

16. S ERIES

16.1. Taylor Series. The Taylor series of a function f (x) at x = a is of the form
∞ f (n) (a)
(x − a)n .
X
f (x) =
n=0 n!
In particular, if a = 0, it is called the Maclaurin series of f (x):
∞ f (n) (0)
xn.
X
f (x) =
n=0 n!

MATLAB has the command taylor to produce the Taylor series of a function up to certain
order:
>> taylor(function, variable, point, 'Order', number)
x2 + 3
For example, find the Taylor series of f (x) = at x = 3 of order < 6:
x +5
>> taylor((x2+3)/(x+5), x, 3, 'Order', 6)
ans = (9*x)/16 + (7*(x - 3)2)/128 - (7*(x - 3)3)/1024 + (7*(x - 3)4)/8192
- (7*(x - 3)5)/65536 - 3/16

16.2. Fourier Series. The fourier series of a periodic function f (x) with periodic 2L is given by
a0 X∞ nπx X ∞ nπx
f (x) = + a n cos + b n sin ,
2 n=1 L n=1 L
where
2 L πnx 2 L πnx
Z Z
an = f (x) cos dx and b n = f (x) sin d x.
L −L L L −L L
There is no simple command in MATLAB to produce the coefficients for fourier series. But
we can define by ourself.
(
0 if − 2 ≤ x < 0,
For example, let f (x) = and f (x) = f (x + 4) for all x; so L = 2.
x if 0 ≤ x < 2,
>> syms x
>> f = piecewise(-2<=x<0, 0, 0<=x<2, x);
>> L = 2;
>> a = @(n) 2/L*int(f*cos(pi*n*x/L), x, [-L, L]);
>> b = @(n) 2/L*int(f*sin(pi*n*x/L), x, [-L, L]);
 INTRODUCTION TO MATLAB 27

We may find a 0 , . . . , a 6 and b 1 , . . . , b 6 using simple commands:

>> a(0:6)
ans = [ 2, -8/pi2, 0, -8/(9*pi2), 0, -8/(25*pi2), 0]
>> b(1:6)
ans = [ 0, 4/pi, -2/pi, 4/(3*pi), -1/pi, 4/(5*pi), -2/(3*pi)]

17. O PERATIONS ON M ATRICES

17.1. Basic Operations. The entries of a matrix shall be entered row by row, while the entries in
each row are separated by a space and the rows are separated by a semi-colon ; . For example,
à !
1 2 3
A= .
4 5 6

>> A = [1 2 3; 4 5 6]
A = 1 2 3
A = 4 5 6
The matrix addition, subtraction and multiplication with scalar can be evaluated using + ,
- and * respectively. For example,
à ! à !
1 2 4 1
A= and B = .
3 4 2 5

>> A = [1 2; 3 4];
>> B = [4 1; 2 5];
>> A + B
ans = 5 3
ans = 5 9
>> A - B
ans = -3 1
ans = 1 -1
>> 3*A
ans = 3 6
ans = 9 12 Ã ! Ã !
4 1 1 2
We illustrate more operations using the previously defined A = and B = .
3 4 2 5
 INTRODUCTION TO MATLAB 28

(i) Transpose AT :
>> A'
ans = 1 3
ans = 2 4
(ii) Rank rank(A):
>> rank(A)
ans = 2
(iii) Determinant det(A), provided that A is a square matrix.
>> det(A)
ans = -2
(iv) Powers An , provided that A is a square matrix. If n < 0, A needs to be invertible (that is,
non-singular).
>> A2
ans = 7 10
ans = 15 22
(v) If A is invertible, its inverse can be evaluated using either A(-1) or inv(A) .
>> inv(A)
ans = -2.0000 1.0000
ans = 1.5000 -0.5000
(vi) The adjoint adj(A), provided that A is a square matrix.
>> adjoint(A)
ans = 4.0000 -2.0000
ans = -3.0000 1.0000
(vii) Matrix product AB , provided that the sizes are matched.
>> A * B
ans = 8 11
ans = 20 23
Moreover, we can generate special matrices using the following command:
(i) Zero matrix 0m×n of size m × n: zeros(m,n) .
>> zeros(2,3)
ans = 0 0 0
 INTRODUCTION TO MATLAB 29

ans = 0 0 0
(ii) Identity matrix In of order n: eye(n) .
>> eye(2)
ans = 1 0
ans = 0 1
(iii) Diagonal matrix with diagonal entries a 1 , . . . , a n : diag([a1 ... an]) .
>> diag(2,3)
ans = 2 0
ans = 0 3

17.2. Row and Column Operations. Let A be a matrix. For example,

 
1 2 3 4 5
 
2 3 4 5 6
A= .
 
3 4 5 6 7
 
4 5 6 7 8

>> A = [1 2 3 4 5; 2 3 4 5 6; 3 4 5 6 7; 4 5 6 7 8];
(i) The size of A is given by
>> size(A)
ans = 4 5
(ii) The (i , j ) entry of A is simply A(i,j) . For example,
>> A(2,5)
ans = 6
(iii) The i th row of A: A(i,:) . For example,
>> A(4,:)
ans = 4 5 6 7 8
If we need more rows, indicate the indices in square brackets. For example,
>> A([2,4], :)
ans = 2 3 4 5 6
ans = 4 5 6 7 8
(iv) The i th column of A is given by A(:,i) .
>> A(:,3)
ans = 3
 INTRODUCTION TO MATLAB 30

ans = 4
ans = 5
ans = 6
(v) If we need more columns, indicate the indices in square brackets.
>> A(:, [3,4])
ans = 3 4
ans = 4 5
ans = 5 6
ans = 6 7
(vi) Submatrix of A. Simply indicate the required rows and columns in square brackets.
>> A([1,2], [3,4])
ans = 3 4
ans = 4 5
We can perform the elementary row operations as follows:
(i) Multiplying the i th row by a constant c: A(i,:) = c*A(i,:) .
>> A(1,:) = -2*A(1,:);
A = -2 -4 -6 -8 -10
A = 2 3 4 5 6
A = 3 4 5 6 7
A = 4 5 6 7 8
(ii) Interchanging the i th and j th rows: A([i,j],:) = A([j,i],:) .
>> A([2,3],:) = A([3,2],:);
A = -2 -4 -6 -8 -10
A = 3 4 5 6 7
A = 2 3 4 5 6
A = 4 5 6 7 8
(iii) Adding c times of the i th row to the j th row: A(j,:) = A(j,:) + c*A(i,:) .
>> A(4,:) = A(4,:) + 2*A(1,:)
A = -2 -4 -6 -8 -10
A = 3 4 5 6 7
A = 2 3 4 5 6
 INTRODUCTION TO MATLAB 31

A = 40 -3 -6 -9 -12
(iv) The reduced row echelon form of A is given by
>> rref(A)
ans = 1 0 -1 -2 -3
ans = 0 1 2 3 4
ans = 0 0 0 0 0
ans = 0 0 0 0 0
Similarly, we can perform the elementary column operations as follows:
(i) Multiplying the i th column by a constant c: A(:,i) = c*A(:,i) .
(ii) Interchanging the i th and j th columns: A(:,[i,j]) = A(:,[j,i]) .
(iii) Adding c times of the i th column to the j th column: A(:,j) = A(:,j) + c*A(:,i) .

18. L INEAR S YSTEM

18.1. Homogeneous Linear System. Let A be a matrix. The solution set of the homogeneous
linear system Ax = 0 can be obtained by back-substitution from its reduced row echelon form.
 
0 3 −6 6 4
 
For example, let A = 
3 −7 8 −5 8.

3 −9 12 −9 6
>> A = [0 3 -6 6 4; 3 -7 8 -5 8; 3 -9 12 -9 6];
>> rref(A)
ans = 1 0 -2 3 0
ans = 0 1 -2 2 0
ans = 0 0 0 0 1
Alternatively, we note that the solution set of Ax = 0 is precisely the nullspace of A. The
command null(A) finds a basis for the nullspace of A (using column vectors).
>> null(A)
ans = -0.7952 0.1467
ans = -0.3797 0.5633
ans = 0.2256, 0.6983
ans = 0.4155 0.4166
ans = -0.0000 -0.0000
 INTRODUCTION TO MATLAB 32

Let v1 = (−0.7952, 0.3797, 0.2256, 0.4155, −0.0000)T and v2 = (0.1467, 0.5633, 0.6983, 0.4166, −0.0000)T .
Then Ax = 0 has a general solution x = c 1 v1 + c 2 v2 .

18.2. Non-Homogeneous Linear System. Let A be an m×n matrix and b an m×1 vector. Recall
that the linear system Ax = b can be solved as follows:

(i) Find a particular solution to Ax = b (if the system is consistent), say xp .

(ii) Find the general solution to the homogeneous system Ax = 0, say xh .

Then the general solution to Ax = b is x = xp + xh .

One may use the back-substitution method to solve the system from the reduced row eche-
lon form of the augment matrix (A | b).
 
−5
 
For example, let A be the matrix previously defined, and b = 
 9 .

15
>> b = [-5; 9; -15];
>> rref([A b])
ans = 1 0 -2 3 0 -24
ans = 0 1 -2 2 0 -7
ans = 0 0 0 0 1 4
Alternatively, both linsolve(A,b) and A\b gives a particular solution to the system Ax =
b, provided that the system is consistent.
>> linsolve(A,b)
ans = -17.000
ans = 0
ans = 3.5000
ans = 0
ans = 4.0000
So the system Ax = b has a particular solution v = (−17, 0, 3.5, 0, 4)T .
Recall that Ax = 0 has a general solution c 1 v1 + c 2 v2 . Then Ax = b has a general solution
x = v + c 1 v1 + c 2 v2 , where c 1 , c 2 are arbitrary parameters.

18.3. Least Squares Solution. Sometimes the system Ax = b is inconsistent, we prefer to get
the least squares solution, i.e., the solution to AT Ax = AT b. A least squares solution can also
be obtained by linsolve(A,b) or A\b .
 INTRODUCTION TO MATLAB 33
   
4 0 2
   
For example, A = 0 2 and b =  0 
  
.
1 1 11
One checks that the system is inconsistent. But we have least squares solution x = (1, 2)T .
>> A = [4 0; 0 2; 1 1];
>> b = [2; 0; 11];
>> linsolve(A,b)
ans = 1.0000
ans = 2.0000

19. M ORE M ATRIX O PERATIONS

19.1. Orthonormal Basis. Let V be a vector space spanned by vectors u1 , . . . , uk . Using Gram-
Schmidt process, one obtains an orthonormal basis for V . In MATLAB, orth can be used to
generate orthonormal basis for V .
More precisely, orth(A) gives an orthonormal basis for the column space of the matrix A.
For example, let V be the vector space spanned by

v1 = (−10, 2, −6, 16, 2), v2 = (13, 1, 3, −16, 1), v3 = (7, −5, 13, −2, −5), v4 = (11, 3, −3, 5, −7).

We shall first define a matrix A whose columns are v1 , v2 , v3 , v4 . Since it is easier to input row
vectors in MATLAB, we may view each vi as a row vector and take transpose:
 
 T −10 13 7 11
v
 1
 
 2 1 −5 3 
v   
 2
A=  =  −6 3 13 −3 .
 
v3   
   16 −16 −2 5 
v4
 
2 1 −5 −7

>> A = [-10 2 -6 16 2; 13 1 3 -16 1; 7 -5 13 -2 -5; 11 3 -3 5 -7]';

>> orth(A)
ans = -0.6140 -0.5638 0.3459 0.3489
ans = 0.0720 -0.0389 0.4409 0.2483
ans = -0.3406 -0.0986 -0.8114 0.3016
ans = 0.7011 -0.6159 -0.1476 0.3117
ans = 0.1016 0.5399 0.0764 0.7928
The columns of the resulting matrix is thus an orthonormal basis for V .
 INTRODUCTION TO MATLAB 34

19.2. Eigenvalues. Let A be a square matrix of order n. Its characteristic polynomial is det(x In −
A). Alternatively, we can use the command charpoly .
 
4 −1 6
 
For example, let A = 
2 1 6 .

2 −1 8
>> A = [4 -1 6; 2 1 6; 2 -1 8];
We shall first declare the variable of the polynomial, e.g., x:
>> syms x;
Then either of the following commands evaluates the characteristic polynomial of A:
>> det(x*eye(3)-A)
ans = x3 - 13*x2 + 40*x - 36
>> charpoly(A,x)
ans = x3 - 13*x2 + 40*x - 36
The roots of the characteristic equation of A are precisely all the eigenvalues of A:
>> solve(charpoly(A,x) == 0, x)
ans = 2
ans = 2
ans = 9
Alternatively, MATLAB can evaluate the eigenvalues using eig .
>> eig(A)
ans = 9.0000
ans = 2.0000
ans = 2.0000
The command [matrix1, matrix2] = eig(A) can return the corresponding eigenvec-
tors. Here matrix2 is a diagonal matrix whose diagonal entries are the eigenvalues of A
(counting multiplicities) and matrix1 is a matrix whose columns are the corresponding (unit)
eigenvectors of A.
>> [P, D] = eig(A)
P = -0.5774 -0.6122 0.3205
P = -0.5774 -0.7873 -0.9112
P = -0.5774 0.0728 -0.2587
D = 9.0000 0 0
 INTRODUCTION TO MATLAB 35

D = 0 2.0000 0
D = 0 0 2.0000
Note that A is diagonalizable if and only if matrix2 is invertible.
In this problem, one verifies that P is invertible and P −1 AP = D :
>> inv(P)*A*P
ans = 9.0000 -0.0000 -0.0000
ans = 0.0000 2.0000 0
ans = 0.0000 -0.0000 2.0000

19.3. Differential Equation in Matrix Form. A linear differential equation may be written in
matrix form. For example,
dx


 = x + 2y + 1,
dt

 dy
= −x + y + t .


dt
It can be expressed in matrix form:
à ! à !à ! à !
x0 3 −4 x 1
= + .
y0 4 −7 y t
à ! à ! à !
x 1 2 1
Let x = ,A= and b = . Then the differential equation is simply x0 = Ax + b.
y −1 1 t
The command dsolve can also be used for differential equation in matrix form:
>> syms x(t) y(t);
>> X = [x; y];
>> A = [3 -4; 4 -7];
>> b = [1; t];
>> [Sx, Sy] = dsolve(diff(X,t) == A*X + b)
Sx = (exp(-5*t)*(C1 + (2*exp(5*t)*(10*t - 7))/75))/2 + 2*exp(t)*(C2 + (exp(-t)*(t
- 1))/3)
Sy = exp(-5*t)*(C1 + (2*exp(5*t)*(10*t - 7))/75) + exp(t)*(C2 + (exp(-t)*(t
- 1))/3)