Ville-Pekka Eronen, Jarkko Peltomäki
Department of Mathematics and Statistics, University of Turku
Last updated: 31.8.2017

This document contains an introduction to the MATLAB software (later simply Matlab). The aim of the material is to give sufficient knowledge on the software to attend to Tero Aittokallio's course Scientific Computing; it is not meant to be an introduction to programming or problem solving with Matlab. The material is mainly for self-study, experimenting is arguably one of the best ways to learn specific software and programming.


Short History of Matlab

Matlab (stands for matrix laboratory) is a proprietary software that is meant for matrix-based numerical computing. It was developed in the 1970's by Cleve Moler who wanted his students too access Fortran-based linear algebra libraries without knowing Fortran itself. Later Moler commercialized the software and cofounded the MathWorks software company in 1984. The company actively develops Matlab even today: the 50th release version, dubbed R2017a, was released in March 2017. What makes Matlab special is that it is a pure matrix language; even real-valued variables are strictly 1-by-1 matrices. A significant portion of numerical tasks is naturally interpreted to the language of matrices, so Matlab, being based on the optimized LINPACK linear algebra library, is very efficient and popular numerical computation software. Its matrix operations are notably faster than, for instance, Mathematica's. Unlike Mathematica, Matlab cannot natively manipulate symbolic computations, even though the feature is available as an additional toolbox. Matlab is a proprietary software, and its licence is quite expensive. Its is, however, possible for students to get a cheaper licence. It is possible to use free software, such as Scilab or GNU Octave, to learn the principles of Matlab. These programs are not complete Matlab clones, but the basic programming covered in this document is very similar. Matlab is available for the usual platforms: Unix, Mac OS X, and Windows.

Links to Other Materials

The User Interface and the Interpreter

By starting the Matlab software, the following window opens:

In the middle of the above picture, there is the interpreter (command window) where you can interactively run Matlab commands. Displayed on the left, there is the current folder. On the right, there is the workspace that records the current variables, i.e., the state. Below it, we see the command history.

The interpreter's prompt >> is distinct from the printed lines not displaying it. A command is executed by pressing the ENTER button. For example:

x = 5 + 1
x =


phi = (1+sqrt(5))/2
phi =


We observe that Matlab prints output after every command. This interactivity is often unnecessary, so printing can be suppressed by ending a command with a semicolon ;. Using a semicolon or a comma, several commands can be written on a single line, if desired. The percent sign % starts a comment. The syntax of Matlab is straightforward to anyone who has used some C-like programming language. We will not dwell on different aspects of the syntax here, it should become evident from the examples and the experiments by the reader. However, we will say this: Matlab is a dynamically typed language. Notice that by creating a new variable, it is possible to overwrite a built-in Matlab command. Defining a variable named sin prevents the usage of the function sin; this should be avoided.

More information on a Matlab command can be obtained by writing help command or by writing the command into the search box (top right in the above picture). By simply writing help, a list of available commands organized by topic is printed. The command lookfor word finds all commands that contain the given word on the first line of their help message (the first line contains a description for the whole function). In addition, the command lookfor word -all lists all commands whose whole help message contains the given word. The web page of Mathworks also provides a wealth of information.

Exercise. Make yourself comfortable with Matlab's user interface and try basic arithmetic operations in the interpreter. Use the help command to study the command format, and print the golden ratio with greater accuracy. Does the amount of displayed digits affect the precision used in the computations? Find out what format compact does.

By writing the above example in the interpreter, we see that the defined variables x and phi will stay in memory and are displayed in the workspace. The command who lists all defined variables, clear deletes all of them, and clear x removes the variable x. Try also the command whos. Like in a calculator, Matlab offers the special variable ans that contains the result of the previous command. The clc clears the text displayed in the interpreter.

The command disp is used for printing. For example:


Printing strings works as follows (in Matlab a string is indicated by single quotes):

disp('Hello world!')
Hello world!

In Matlab, strings are arrays of characters, so more complex output can be generated as follows (cf. the creation of vectors and matrices below).

disp(['The value of the variable phi is ' num2str(phi) '.'])
The value of the variable phi is 1.618.

The elementary functions, such as abs, sin, cos, tan, exp, log, sinh, are built in Matlab. More information can be obtained by typing help elfun. Special functions, like Bessel functions, are found by entering help specfun.

Exercise. Visit the Matlab documentation on page Complex Numbers and find out how complex numbers work in Matlab. Try calling elementary functions with complex parameters.


Like mentioned already, Matlab is a matrix language. The easiest way to define a matrix is to give its elements by rows:

A = [1 2 3; 4 5 6; 7 8 9]
A =

     1     2     3
     4     5     6
     7     8     9

The definition can be split onto several rows. Then the semicolon is not needed.

A = [1 2 3
4 5 6]
A =

     1     2     3
     4     5     6

By writing three dots, a long row can be split onto two separate lines:

A = [1 2 3 4 5 6 ...
7 8 9
10 11 12 13 14 15 16 17 18]
A =

     1     2     3     4     5     6     7     8     9
    10    11    12    13    14    15    16    17    18

The elements of a matrix can be separated by comma instead of a space. Naturally the elements can be defined as results of arithmetic operations.

A = [1, 2*sin(pi/4), exp(2)]
A =

    1.0000    1.4142    7.3891

The command length returns the length of a vector:

ans =


On the other hand, the dimension of a matrix is given by the command size, try it! Matlab has several built-in commands for matrix generation. For example, the command eye returns an identity matrix:

ans =

     1     0     0
     0     1     0
     0     0     1

Exercise. Use the help command, and get to know the commands rand, ones, and zeros.

The usual arithmetic operations +, -, *, /, \, ^, ', addition, subtraction, multiplication, division from right and left, exponentiation, transposition, are available for matrices. The programmer needs to make sure that the operated matrices are of the appropriate size. It is easier to get to know the operations by experimenting; we shall return to the division later when we discuss solving a system of equations. It is worth mentioning here that adding or subtracting a scalar has special syntax:

A = eye(3)
A - 1
A =

     1     0     0
     0     1     0
     0     0     1

ans =

     0    -1    -1
    -1     0    -1
    -1    -1     0

We see that Matlab can automatically interpret the subtracted scalar to be multiplied by the matrix ones(3). It is also worth mentioning, that ' is Hermitian transpose, that is, in addition to the usual transposition, also complex conjugation is applied. The usual transposition is available with the command transpose or with the notation .'.

Multiplication, exponentiation, and division also have elementwise versions. So for example, .* computes the elementwise product of two matrices instead of the matrix product. For example:

A = [1 2; 3 4], B = [1 1; 2 2]
A .* B
A =

     1     2
     3     4

B =

     1     1
     2     2

ans =

     1     2
     6     8

Many Matlab functions, such as the elementary functions mentioned above, operate elementwise. For instance, the function sin computes the sine of every element of a matrix:

ans =

    0.8415         0
         0    0.8415

Colon Notation and Indexing

It is often necessary to produce vectors where the distance between two consecutive elements is constant. This is easily accomplished with the colon notation:


This produces a vector containing the numbers start, start + step, start + 2* step, ... that are smaller than or equal to end. For example:

ans =

  Columns 1 through 7

    1.0000    1.1000    1.2000    1.3000    1.4000    1.5000    1.6000

  Columns 8 through 11

    1.7000    1.8000    1.9000    2.0000

ans =

    1.0000    1.3000    1.6000    1.9000

If step is omitted, then its value is assumed to be 1. The following code produces the same result using the for-loop. This method should not be used, this is only for introducting the syntax of the for-loop and for showing how to add elements to vectors.

X = [];
s = 1; e = 2; t = 0.1;
for n = 0:(e-s)/t
    X = [X s + n*t];
X =

  Columns 1 through 7

    1.0000    1.1000    1.2000    1.3000    1.4000    1.5000    1.6000

  Columns 8 through 11

    1.7000    1.8000    1.9000    2.0000

Exercise. Figure out how the command linspace works and how it is different from the colon notation. The command logspace might also be useful later.

The next code shows how to select specific elements, rows, and columns from a matrix.

A = [1 2 3 4; 5 6 7 8; 9 10 11 12]
A(2,3)                                   % The element at (2,3).
A(8)                                     % The element at (2,3)
A(:,1)                                   % 1st column
A(1,:)                                   % 1st row
A(:)                                     % The whole matrix as a column vector
A(2:3,2:3)                               % Submatrix, rows 2-3, columns 2-3
A(:,[1,4])                               % All rows, columns 1 and 4.
B = logical([0 1 0 0; 1 0 1 0; 0 0 1 0]) % Creates a "logical" matrix which tells which elements of A are to be chosen.
A(B)                                     % Selects the elements of A whose corresponding positions in B contain 1, ordering is top to bottow, left to right.
A =

     1     2     3     4
     5     6     7     8
     9    10    11    12

ans =


ans =


ans =


ans =

     1     2     3     4

ans =


ans =

     6     7
    10    11

ans =

     1     4
     5     8
     9    12

B =

     0     1     0     0
     1     0     1     0
     0     0     1     0

ans =


A matrix can be analogously modified:

A(2,2:3) = [0 0]    % Change the columns 2 and 3 on row 2.
A(:,4) = [12; 8; 4] % Change 4th column.
A =

     1     2     3     4
     5     0     0     8
     9    10    11    12

A =

     1     2     3    12
     5     0     0     8
     9    10    11     4

Exercise. Define some matrix and add its first row multiplied by 3 to its second row. The colon notation thus allows to implement the elementary row operations for matrices. It is thus quite straightforward to program a program that puts a matrix into row echelon form.

On Strings

Strings can be searched and edited with regular expressions. The next example finds the nonoverlapping occurrences of the string aa in a longer string.

regexp('abbabaabaaa', 'aa')
ans =

     6     9

Naturally it is possible to write proper regular expressions. The next example finds "greedily" all occurrences of words that begin with k and end with o from a longer string.

s = 'koko kokoa koko kokko kokoon';
regexp(s, 'k[\w]*o')                 % returns the starting positions of the matching occurrences
regexp(s, 'k[\w]*o', 'match')        % returns the matching words themselves
ans =

     1     6    12    17    23

ans = 

    'koko'    'koko'    'koko'    'kokko'    'kokoo'

Substrings are replaced as follows:

regexprep(s, 'k', 'o')    % change letters k to letters o
ans =

oooo ooooa oooo ooooo ooooon

Read more on regular expressions on the page Match regular expression (case sensitive).

m Files and On Writing Own Functions


Often there is a need to save the typed commands to a file after an interactive session. This is accomplished by saving the code into a text file with extension .m. The code can be written with any editor, but it is very convenient to do it using Matlab's built-in editor. The editor can be opened by selecting New -> Script from the user interface. The same can be accomplished by typing edit filename into the interpreter. Let us copy the following code to the editor and save it as the file testi.m into the current folder (so that Matlab finds it, this is important).

x = randn(1, 100);
m = mean(x)
s = std(x)

After this, write testi in the prompt (you can also click Run in the editor).

m =


s =


The result is thus the same as if the contents of the file testi.m were written in the interpreter. Notice that exactly the same output was produced and after the execution the defined variables x, m, and s stayed in the workspace. By using m files, it is thus possible to initialize variables with proper values in the beginning of an interactive session, for instance. It is possible in an m file to call the code of another m file in the exact same way. Remember that the contents of a m file can be commented using the comment sign %; this sould always be done. If you select some lines of code in the editor, their contents can be executed by pressing the F9 button.

An m file prepared in this fashion is called a script, and it cannot take any parameters when called. The next subsection tells that m files may contain functions, which naturally can take parameters.

Matlab can be used to produce documents that contain code, text, and other markup and that can be compiled into a HTML or PDF file. This document itself is a single m file that is compiled into a HTML file using the publish command. More information on this on the page Publishing Markup.

Function m Files

Function definitions are written to m files with the following syntax:

function y = avg(x)
y = sum(x)/length(x);

In order for Matlab to recognize the m file as a function definition, the keyword function needs to be found on the first line of the m file. The m file containing the function must have the same name as the function (here avg.m) so that Matlab finds the desired function. Notice that the commands within a function should end with a semicolon; otherwise lots of unnecessary output might be produced. Remember that, as is usual with common programming languages, variable definitions inside a function definition are local and are not available after the execution of the function.

In the above function, y is the value returned by the function, and it is set on the file line of the function definition. The name of the function is avg and its parameter's name is x. A function may take several parameters: avg(x,y) etc.

A function may also return several values. Then the syntax is the following:

function [a,b] = foobar(x)
a = sum(x);
b = sum(x)/length(x);

This function is called as follows:

[a,b] = foobar(x)

If the user does not care for all the values, then only the first one can be saved:

a = foobar(x)

The preceding bracket notation can be analogously be used to pick one or more return values of a function. If a function returns nothing, the return value can be omitted:

function func(x)

Matlab supports optional parameters, that is, not all parameters need to be specified. Then it is up to the function if this raises an error or if it modifies the behavior of the function in some way (cf. More on Programming Functions below).

A single m file may contain the definition of two or more functions, but these extra functions are local to the m file. Only the file's "main" function may call these auxiliary functions. In this case, all function definitions must be terminated with the keyword end; otherwise it would be unclear where each function definition ends.

Anonymous Functions

Several Matlab functions take functions as parameters, like the function fzero that finds a root of a function. Let us suppose that the file fun.m has the following contents:

function y = fun(x)
y = 1./(exp(x) + sin(x)) - x.^2;

The function defined in the file thus computes the value of the function $\frac{1}{e^x + \sin{x}} - x^2$. Let us find its root in the interval $[-0.5,1]$:

fzero(@fun, [-0.5, 1])
ans =


The significance of the @ sign is that it creates a reference (function handle) to the function fun. References can be saved into variables, so the following code works as expected:

f = @fun;
fzero(f, [-0.5, 1])
ans =


It is tedious to write a separate m file for each little function, so Matlab provides syntax for defining a so-called anonymous function. The following code does the same as the preceding code without the need to write a separate m file:

f = @(x) 1./(exp(x) + sin(x)) - x.^2;
fzero(f, [-0.5, 1])
ans =


An anonymous function can also depend on a parameter:

a = pi;
f = @(x) 1./(exp(x) + sin(x)) - a*x.^2;
fzero(f, [-0.5, 1])
ans =


Read more on anonymous functions here.

The Location of m Files

By default Matlab searches for scripts and functions from the directories listed by the path command. The first directory in the list is the current working directory. The function addpath can be used to add more directories to the list. The directories matter especially when several m files with the same name are specified. In this case, Matlab picks the first m file it encounters when searching the directories given by addpath in order. The command which tells which m file Matlab attempts to use. Note that an m file should not have the same file as a built-in Matlab command as the builtins are always called first.

Matlab Toolboxes

Matlab provides several toolboxes. These are typically collections of functions, that is, collections of m files, but they can also contain applications and data.

Toolboxes available during a session are found with the command ver. By executing the command, we learn, for instance, that the "optimization toolbox" is available. Thus all functions provided by the toolbox are directly available. This page lists additional purchasable toolboxes. Matlab Central provides third-party toolboxes; consult their readme files to see how to install them. If a toolbox is simply a directory containing m files, it is sufficient to add it using the addpath command. If the directory G:\ELBM contains m files, they are made available as follows:

help ELBM          % Matlab does not know this function.
help ELBM          % Now it does.

More on Programming Functions

Matlab naturally has the usual control structures if, for, while, etc. We will not look at these in this document; for more, see the page Control Flow. In some of control structures it is useful to use logical operations. These are found with command help ops. The implementation of several Matlab commands are viewable using the type command. Try type('quad').

Next we introduce a function that computes an approximation of pi based on the series $4(\frac11-\frac13+\frac15-\frac17+\cdots)$ (very slow convergence). The behavior of this function depends on the number of parameters given; it also prints an error if a parameter is invalid.

function p = pi_approx(iterations)
  % If the user did not specify the number of iterations, do 1000 iterations.
  if nargin() ~= 1
    iterations = 1000;
  if iterations <= 0
    error('The number of iterations must be a positive integer.')

  s = 0;
  for n = 0:(iterations-1)
    if mod(n, 2) == 0
      s = s + 1/(2*n + 1);
      s = s - 1/(2*n + 1);

  p = 4*s

The command nargout returns the number of return values requested by the caller. This is sometimes good to know if some return value is expensive to compute; if it is not requested, it is useless to compute it.

Acquiring input from the user is performed with the command input. This command takes as a parameter a string to be displayed when the user is prompted; it returns the value inputted by the user. The function below is the above function modified; now it asks the user to input the number of iterations.

function p = pi_approx_inter()
  iterations=input('How many iterations? ');
  if ~isnumeric(iterations)
       error('The number of iterations must be numeric')
  if iterations <= 0
    error('The number of iterations must be a positive integer.')

  s = 0;
  for n = 0:(iterations-1)
    if mod(n, 2) == 0
      s = s + 1/(2*n + 1);
      s = s - 1/(2*n + 1);

  p = 4*s

As Matlab is dynamically typed, the functions shold be programmed to accept different types of parameters. For instance, a function computing the value of a real function at a given point might as well accept a vector or matrix of points as its input. Accomplishing this in Matlab is often straightforward due to the elementwise operations on matrices. The next function computes the value of the sin function, using the Taylor series, for all elements of the input matrix (or vector, or number).

function y = sin_taylor(A, iterations)
  if nargin() ~= 2
    iterations = 10;
  if iterations <= 0
    error('The number of iterations must be a positive integer.')

  y = 0;
  for n = 0:(iterations-1)
    if mod(n, 2) == 0
      y = y + (1/factorial(2*n+1))*A.^(2*n+1);
      y = y - (1/factorial(2*n+1))*A.^(2*n+1);

With the default number of iterations, the result agrees with the output of Matlab's sin function.

sin_taylor([0 0.1; 0.5 1]), sin([0 0.1; 0.5 1])
ans =

         0    0.0998
    0.4794    0.8415

ans =

         0    0.0998
    0.4794    0.8415

Exercise. Program a function that takes as its input a square matrix A (check that the user does not give an inappropriate matrix) and computes the matrix $e^A$ with the series $e^A = A + 1/2!A^2 + 1/3!A^3 + \ldots$ using as many terms as the user desires.

On Debugging

If the code contains an error, Matlab prints a reason for it. Sometimes Matlab prints the line number where the error occured. Click the link provided by Matlab to get to edit the errorneous line. An error message can also be studied with the command lasterror.

err = lasterror

The code can be debugged and single-stepped by using Matlab's debugger. The debugger is activated by adding a breakpoint into the code. This can be done either by clicking the dash next to a line number (it turns into a red ball) or by the dbstop command. Let us try this in the case of the following code from Matti Vuorinen's course (lecture notes p. 64):

dbstop in e310.m at 15

This sets a breakpoint at line 15. When executed, Matlab goes into debugging mode, which is recognizable from the letter K next to the prompt. Matlab commands can be used normally in debug mode, and the code can be single-stepped. The next line is executed with the dbstep command. The command dbcont resumes the execution until the next breakpoint. The debug mode can be exited with the command dbquit; this terminates the execution of the code.


By using the debugger, it is thus possible to check that the variables have the values they should have.

Sometimes valid code executes slowly. With the commands tic and toc you can measure how long the code spends time in a section of code.

for i=1:10000
Elapsed time is 0.001604 seconds.
Elapsed time is 0.000280 seconds.

More time-related functions are found by typing help timefun.

Importing Data to Matlab

Matlab has many commands for importing date such as load, fscanf, textscan, importdata, and xlsread. More functions are found from the page Supported File Formats for Import and Export.

Let us first look at the command load. It can load mat files and ASCII files containing numerical values. Mat is Matlab's own format, which is used to save workspace variables with the save command.

m = 1;
M = [1,2; 3,4];
save('test.mat', 'M');  % Save matrix M.
who;                    % Look which variables occupy workspace.
who;                    % Now the workspace is empty.
who;                    % M loaded.

If you do not specify the variables to be saved, the save command saves the state of the whole workspace. If no filename is specified, the default name matlab.mat will be used.

The file plotdata.txt contains two columns of numerical values. These can be loaded as a matrix by typing

M = load('plotdata.txt');

Occasionally the data does not consist of only numerical values. Then other functions can be used to import the data. Let us next look at importing xlsx data with the xlsread command.

xlsx file

The data consists of three sheets. In the first sheet there is a simple table. Let us first check how we can read this by typing help xlsread in Matlab.

Read the table with the following command

M = xlsread('xlsxtest.xlsx')

We observe that no strings are read and the numerical data is stored in the variable M. The strings can be saved to the variable R as follows:

[M,R,U] = xlsread('xlsxtest.xlsx')

The variable U contains "raw data". Its class is cell array and each cell includes one read element. In general a cell can contain any data type. The data contained in a cell can be extracted by using braces (look the exqample below) or by using the function cell2mat.


Let us move on to the next sheet. In this data there are two missing values and one NA value. In addition, there is a new column of strings. Let see how xlsread will deal with this.

[M,R] = xlsread('xlsxtest.xlsx',2)

We observe that the missing values and the NA value has been replaced by NaN (double), which stands for "not a number". Furthermore, the string column was correctly read in the variable R.

In the third sheet there are decimal commas. These are nontrivial in the sense that Matlab uses points to indicate decimals.

[M,R] = xlsread('xlsxtest.xlsx',3)

We observe that xlsread was able to change decimal commas to dots correctly.

Exercise: Try to read the following datas by using importdata function. (spoiler:problems ahead)

dotdata specialdata commadata

The import wizard (File->Import Data...) of Matlab uses the function importdata. As can be seen from the exercise it is not always straightforward to use and sometimes it is necessary to use other functions.

One easy way to import the datas in the exercise is to read them in excel, save an xlsx-file and use the function xlsread. Another way is to use the function textscan or fscanf. For example, the data in the file "dotdata" can be read as follows:

fo=fopen('G:\matlab\dotdata.txt'); % open the file.
header=fscanf(fo,'%s %s %s %s',4); % variable for the header row.
m=textscan(fo,'%s %f %f %f %f');   % The rest of the rows are of the same form.
fclose(fo);                        % close the file.
M=[m{2},m{3},m{4},m{5}];           % Creating the numeric matrix.

We can also change the decimal commas into decimal dots in the file "commadata.txt" with fscanf.

commadata=fscanf(fo,'%c') % '%c' preserves empty spaces
dotdata=regexprep(commadata,',','.') % Replace commas by dots in the string pilkkudata.
fc=fopen('G:\matlab\commadatamod.txt','w');  % open and, if necessary, create the file "pilkkudatamuok".

More help for importing data can be found on page

Exporting Data from Matlab

For most of the commands to import data there are corresponding commands to export data. An example of this correspondence are the save and write commands. Detailed information on commands to export data is found from the two preceding links. Let us as an example save the matrix

M = [1,2,3; 4,5,6; 7,8,9]

as an excel file. The right syntax is found with the command help xlswrite.

xlswrite('xlswritetest.xlsx', M);

Exercise. Using the command fprintf, save the above matrix M into a text file in such a way that the columns have headers "one", "two", and "three".

2D Graphics


The command plot draws the given points on the plane and connects the points by a line (by default). For example:

x = [0 1 2];
y = [1 -1 0.5];
plot(x, y);

Matlab sets the x-axis and y-axis limits based on the plotted data. This can be adjusted as follows:

axis([-1,1,-2,2]) % x-axis -1 - 1, y-axis -2 - 2

The scale of the axes is controlled with the commands axis equal (same scale on all axis) and axis square (axes are equally long on the screen). Resetting is done by typing axis normal. Notice that these commands need to be given after the plot command (Matlab cannot know the data before it has drawn it). The command grid displays a grid on the plot. Try it!

The commands axis and grid implicitly modify the current figure. Matlab may have several figures open at the same time. A new figure is opened with the command figure; they are referred to by a running sequence of integers. The command figure(n) changes the current figure, which means that the commands axis, grid, etc. operate on another figure. In the next example, we draw the sine and cosine functions in different figures with different scales.

x = linspace(0, 2*pi);
y = sin(x); z = cos(x);
figure;                  % Create new figure, we are now modifying implicitly figure no. 1.
plot(x, y);
axis equal;
figure;                  % Modifying implicitly figure no. 2.
plot(x, z, 'x');         % Draw only data points, not connecting lines.
axis([0, 2*pi, -2, 2]);
figure(1);               % Modifying again figure no. 1.

The figure can be saved in common formats with the command save. For instance, print('image', '-dpdf') saves the current figure as a PDF file (saved in the working directory as the file image.pdf). The parameter -deps saves the figure as EPS. Saving is also accomplished by selecting File -> Save As from the menu of the figure window.

It is possible to draw several plots in the same figure. Calling the command plot erases the previous picture. To prevent this, use the command hold on. Now all calls to plot draw on top of what was previously displayed until the command hold off is used. The plot command can handle several data sets at once:

  plot(x1, y1, 'string1', x2, y2, 'string2', x3, y3, 'string3', ...)

The given strings (LineSpec) control the line color and shape. Try the values '--or' and ':^g'. More information on this on the page Plot.

Exercise. Give a vector of complex numbers for the command plot. Using the above link, find out how to draw the data points without the connecting line segments.


The preceding plot command drew the data points given by two vectors and connected them by line segments. Indeed, this command merely draws the given data, which does not need to represent a function. If you want to plot a built-in Matlab function or a function you have written yourself, the command fplot can be used instead.

fplot(@(x) x.^2 .* sin(x), [0, 4*pi]);

Now Matlab automatically handles the discretization of the interval.

Plotting parametrized functions also works:

fplot(@(t) cos(3*t), @(t) sin(2*t));

Naturally, the usage of polar coordinates is possible:

t = linspace(0,2*pi);
polar(t, sin(2*t).^2);

Multiple Plots in the Same Figure

The command subplot divides the figure into multiple subplots. Calling subplot(m, n, p) divides the figure into m*n subplots that are arranged as a m*n grid. The parameter selects the pth subplot as the current figure (numbering is by rows). In the next example, we create two horizontally adjacent plots, draw to them, and also set titles for the subplots (study the command title in more details).

close all;                    % Close all figure windows.
subplot(1, 2, 1);
fplot(@(x) sin(x), [0, 2*pi]);
title('Sine Curve');
subplot(1, 2, 2);
fplot(@erf, [-1,1]);
title('Something Else');


Here is one more example on adding text to plots, adjusting line width etc. Try it out!

x = linspace(-2*pi, 2*pi); y1 = sin(x); y2 = cos(x);
plot(x, y1, x, y2)
legend('y = sin(x)', 'y = cos(x)');
xlabel('text below the x-axis alle');
ylabel('text next to the y-axis');
axis([-5, 5, -1.5, 1.5]);
set(gca, 'FontSize', 15);           % gca refers to the axes of the current figure
set(gca, 'FontWeight', 'bold');

3D Graphics

The command plot3 is analogous to the command plot:

x = [-1:0.5:1]; y = [2 2 2 3 3]; z = [1:5];
plot3(x, y, z);
xlabel('x'); ylabel('y'); zlabel('z');

Plotting a space curve also works in the same way:

fplot3(@(t) sin(t), @(t) cos(t), @(t) t, [0, 8*pi]);

Let us next consider drawing a surface (that is, a plot of a function from the plane to the line). This is accomplished with the surf command. Let us first look at the command meshgrid. This command takes as its parameters two vectors x and y and forms two matrices X and Y. These matrices correspond to the lattice points defined by the vectors x and y.

x = [0:4]; y = x;
[X, Y] = meshgrid(x, y)
X =

     0     1     2     3     4
     0     1     2     3     4
     0     1     2     3     4
     0     1     2     3     4
     0     1     2     3     4

Y =

     0     0     0     0     0
     1     1     1     1     1
     2     2     2     2     2
     3     3     3     3     3
     4     4     4     4     4

Now the elements X(i,j) and Y(i,j) at (i,j) form the lattice point determined by x(j) and y(i). It is easy to evaluate a function in these lattice points. When we provide the parameters X, Y, and the evaluated values to the command surf, Matlab draws the plane pieces connected by the lattice points. Let us as an example draw the plot of the function $\sin{\sqrt{x^2+y^2}}/\sqrt{x^2+y^2}$.

x = linspace(-10,10); y = x; % the limits of x-axis and y-axis are the same
[X, Y] = meshgrid(x, y);
d = sqrt(X.^2 + Y.^2);
Z = sin(d)./d;
surf(X, Y, Z);

Exercise. Draw the preceding function (or your personal favorite function). Open the Tools menu of the figure window, and try out especially the options Rotate 3D and Data Cursor.

Exercise. Use the preceding values, and call contour(X, Y, Z). What is the meaning of the displayed plot? Consider also the command surfc.

Other Types of Graphics

Here we provide a cursory glance on how to visualize data as histograms, pie diagrams etc. See the page Pie Charts, Bar Plots, and Histograms for completeness.

Let us generate a thousand normally distributed numbers, and draw a histogram.

x = randn(1000, 1);

Let us draw a pie diagram out of the numbers given by the vector x (the areas rotate counter clockwise starting at noon).

x = [0.3 5 2.3 0.88 1];

Linear Algebra

Let us consider solving a linear system of equations. Here we will be using Matlab's operations / and \, that is, division from the right and the left. The system $Ax = y$ is solved by multiplying from the left with the inverse of A. In Matlab, this is accomplished by "dividing" by A from the left:

A = [1 2 3; 0 -4 1; 0 3 -1];
y = [1 2 3]';                % column vector
x = A\y
A*x                          % check for validity
x =


ans =


Naturally, right division solves the system $xA = y$. This division operation in Matlab is not the same as multiplying with the inverse, but it has to do with solving the system numerically, cf. mldivide and mrdivide. The division operation gives an answer ever in the case that the system is overdetermined:

A = [1 2 3; 0 -4 1; 0 3 -1; 0 2 1];
y = [1 2 3 4]';
x = A\y
x =


Here Matlab finds the solution in the sense of least squares.

The coefficient matrix should preferably be invertible. Let us have look what happens if it is not.

A = [1 -2 4; 0 0 6; 0 0 1];
det(A)                      % the determinant of A is 0
y = [1 2 3]';
ans =


Warning: Matrix is singular to working precision. 

ans =


As we can see from the output, Matlab gives a nonsensible result, which hints that the system does not have a solution. If the system had a solution, it would belong to the column space of the matrix A. By adding the vector y as a column of the matrix A, we see that the rank increases:

rank([A y])
ans =


ans =


Thus y does not belong to the column space of A, so the system is contradictory. If a system $Ax = y$ has a solution when $A$ is not invertible, then a solution can be found using the Moore-Penrose inverse:

A = [1 3 7; -1 4 4; 1 10 18];
y = [5 2 12]';
x = pinv(A)*y
x =


ans =


If an exact solution does not exist, the above method yields a solution in the sense of least squares. For more on solving systems of linear equations, see here.

The inverse of a matrix is given by the command inv:

A = [1 1; 0 2];
B = [1 0; 5 2];
B/A             % this also works
ans*A           % check for validity
ans =

    1.0000   -0.5000
    5.0000   -1.5000

ans =

    1.0000   -0.5000
    5.0000   -1.5000

ans =

     1     0
     5     2

Above we introduced the commands det and rank that compute the determinant and the rank of a matrix. Eigenvalues and eigenvectors are found with the command eig:

A = [0 -6 -1; 6 2 -16; -5 20 -10];
eig(A)                             % find only eigenvalues (faster)
[V, D] = eig(A)                    % finds eigenvectors (the columns of the result) ja eigenvalues as a diagonal matrix
V*D*inv(V)                         % this gives back the matrix A (diagonal representation)
ans =

  -3.0710 + 0.0000i
  -2.4645 +17.6008i
  -2.4645 -17.6008i

V =

  -0.8326 + 0.0000i   0.2003 - 0.1394i   0.2003 + 0.1394i
  -0.3553 + 0.0000i  -0.2110 - 0.6447i  -0.2110 + 0.6447i
  -0.4248 + 0.0000i  -0.6930 + 0.0000i  -0.6930 + 0.0000i

D =

  -3.0710 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i  -2.4645 +17.6008i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i  -2.4645 -17.6008i

ans =

   0.0000 + 0.0000i  -6.0000 + 0.0000i  -1.0000 + 0.0000i
   6.0000 - 0.0000i   2.0000 + 0.0000i -16.0000 + 0.0000i
  -5.0000 - 0.0000i  20.0000 + 0.0000i -10.0000 + 0.0000i

If the matrix is non diagonalizable, its Scur decomposition (unitary conjugation to an upper triangular matrix) is obtained with the command schur:

A = [13 8 8; -1 7 -2; -1 -2 7];
[U, T] = schur(A)
U =

    0.9428   -0.3333         0
   -0.2357   -0.6667   -0.7071
   -0.2357   -0.6667    0.7071

T =

    9.0000  -12.7279   -0.0000
         0    9.0000   -0.0000
         0         0    9.0000

Other typical decompositions can be found on the page Matrix Decomposition.

Numerical Integration and Differentiation

The command diff returns the differences of the elements of a vector:

x = [4 3 0 -1 5];
ans =

    -1    -3    -1     6

Thus this command can be used to numerically differentiate (naively):

step = 0.1;                      % step size
x = 0:step:5;
f = x.^2 + 3*x;
df = diff(f)/step;               % difference quatient
plot(x, f, x(1:length(df)), df); % the vector df is one element shorter

It is more straightforward to use the gradient command that can also compute the gradient of a function of several variables.

gradient(f, step);

Numerical integration is done with the integrate command. This command is only available after Matlab version R2012A. With an old version, the trapz command below needs to be used.

integral(@(x) x.^2, 0, 4)
ans =


Exercise. Look up the help on the command integral2 and integrate the function $\sin{\sqrt{x^2+y^2}}/\sqrt{x^2+y^2}$ over the area. $[-1, 1] \times [0, 2 \cdot \pi]$.

The command integral can only integrate a Matlab function, it does not accept a data vector as an input. The integral of an unknown function represented only as a data vector can be computed with the command trapz (trapezoidal rule). Let us apply the command to compute the integral of the function $x^2$ over the interval $[0,4]$.

step = 0.1;
x = 0:step:4;
f = x.^2;     % generate the data
trapz(x, f)
ans =


Solving Differential Equations Numerically

Next we consider how to use Matlab to solve differential equations or systems of differential equations. The command used in this document to solve differential equations is ode45, which is recommended as the default option. Matlab has several other solvers for different uses, cf. Choose and ODE Solver.

Matlab only solves ODE's of the form $y' = f(t, y)$, higher-order equations need to by transformed into this form (cf. lecture notes on differential equations or the above link). The command ode45 takes as its parametrs the function on the right side of the above equation, the integration interval, and an initial value. If we want to find a numerical solution to the ODE $y' = 3ty$ in the time interval $[0,1]$ with initial value $y(0) = 1$, then we need to write

[t, y] = ode45(@(t, y) 3*t*y, [0 1], 1);

The command returns two vectors: the solution y as a function of the time t. The integral curve can thus be immediately drawn:

plot(t, y, '-o')

The general solution to the above ODE is $y = C e^{3t^2/2}$. Thus in this case Matlab's solution looks correct. The command ode45 uses the Runge-Kutta method to solve the ODE.

Exercise. Draw the integral curves of the solutions to the ODE $y' = \frac{A}{B}ty$ with different values for the parameters $A$ and $B$. Use an anonymous function that depends on a parameter.

Let us consider next solving a system of differential equations and use as an example solving the second order ODE $y'' = -sin(y)$. Using the substitution $y = y_1$, $y_2 = y_1'$, we see that we need to solve the system $y_1' = y_2$, $y_2' = -sin(y_1)$. For this system, let us write the following fuction (an anonymous function would also do):

function r = dy_en(t, y)
  % Observe that the parameter y is a vector of two elements because the
  % system has two equations. The command od45 requires that the return
  % value is a column vector, so we cannot simply write
  % 'function [y1, y2] = dy(t, y)' as the function definition.
  r = zeros(2, 1);
  r(1) = y(2);
  r(2) = -sin(y(1));

The solution is obtained analogously:

[t, y] = ode45(@dy_en, [0 5], [0 1]);
plot(t, y(:,1), t, y(:,2));

On Minimizing a Function

Commands suitable for minimizing a function can be found from the optimization toolbox. Its documentation is available here or by typing "optimization toolbox" into the help.

Let us seek a minimum of the function $sin(x) + cos(x)$. At least the commands fminsearch and fminbnd of the toolbox are suitable for the task. Documentation tells us that the syntax for fminsearch is the following:

[x, fval, exitflag, output] = fminsearch(fun, x0, options);

The first parameter is the target function, and the second one is the initial search point. Using the parameter options, it is possible to control how much information is printed.

opt = optimset('Display', 'iter');       % Set options to display informations on iterations.
sincossum = @(x) sin(x)+cos(x);
[x, fval, eflag, output] = fminsearch(sincossum, 0, opt);
 Iteration   Func-count     min f(x)         Procedure
     0            1                1         
     1            2                1         initial simplex
     2            4           0.9995         expand
     3            6         0.998499         expand
     4            8         0.996494         expand
     5           10         0.992472         expand
     6           12          0.98438         expand
     7           14         0.968009         expand
     8           16         0.934527         expand
     9           18         0.864728         expand
    10           20         0.714808         expand
    11           22         0.382525         expand
    12           24        -0.333554         expand
    13           26         -1.34737         expand
    14           28         -1.38509         contract outside
    15           30         -1.41225         contract inside
    16           32         -1.41225         contract inside
    17           34         -1.41412         contract inside
    18           36         -1.41412         contract inside
    19           38          -1.4142         contract inside
    20           40         -1.41421         contract inside
    21           42         -1.41421         contract inside
    22           44         -1.41421         contract inside
    23           46         -1.41421         contract inside
    24           48         -1.41421         contract inside
    25           50         -1.41421         contract inside
    26           52         -1.41421         contract inside
    27           54         -1.41421         contract inside
Optimization terminated:
 the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04 
 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04 

ans =


ans =


In the preceding example, the class of the variables opt and output is struct. A struct is essentialy a named list of elements. Above output.funcCount prints those elements of the variable output that are under the name funcCount. In this particular example, there was only one element, but in general there can be several. In the general case, referring to a single value works as below.

structure = struct('names', {'x','y';'z','w'}, 'values', {1,2;3,4});
ans =


ans =


ans =


ans =


ans =


ans =


The command fminbound is used almost in the same way, this time instead of specifying an initial point, we need to specify an interval the minimum is searched for.

[x, fval, eflag, output] = fminbnd(sincossum, 0,2*pi, opt);
 Func-count     x          f(x)         Procedure
    1        2.39996   -0.0618786        initial
    2        3.88322     -1.41286        golden
    3        4.79993    -0.908745        golden
    4        3.88982     -1.41324        parabolic
    5        3.92941     -1.41421        parabolic
    6          3.927     -1.41421        parabolic
    7        3.92696     -1.41421        parabolic
    8        3.92703     -1.41421        parabolic
Optimization terminated:
 the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04 

ans =


ans =


Finding Roots

Finding a root of a one-variable function is accomplished with the fzero command. Its parameters are a reference to a function and an initial point. This command was earlier introduced in the section on anonymous functions.

fsini = @(x) sin(x);
[x, fval, ef, output] = fzero(fsini, 2)
fsini(x)                                % evaluate and check
x =


fval =


ef =


output = 

    intervaliterations: 10
            iterations: 6
             funcCount: 27
             algorithm: 'bisection, interpolation'
               message: 'Zero found in the interval [0.72, 3.28]'

ans =


Instead of specifying an initial point, an interval containing the root can be specified; see help.

The command fzero seeks only one root. All roots of a polynomial function can be found with the command roots.

a = [1, 2, 3, 4, 5, 6]; % specify the coefficients, x^5 + 2*x^4 + ... + 6
x = roots(a)
polyval(a, x)           % evaluate the polynomial at the found points
x =

   0.5517 + 1.2533i
   0.5517 - 1.2533i
  -1.4918 + 0.0000i
  -0.8058 + 1.2229i
  -0.8058 - 1.2229i

ans =

   1.0e-13 *

   0.3464 + 0.1821i
   0.3464 - 0.1821i
  -0.1332 + 0.0000i
   0.1421 + 0.1998i
   0.1421 - 0.1998i

Exercise. Find some root of the above polynomial using the command fzero.

The command fsolve is used to find a common root of several functions, that is, it can be used to solve a nonlinear system of equations.

nlineq = @(x) [sum(x.^2)-4, [-1,2]*transpose(x.^2)];
opt = optimset('Display', 'iter');
x = fsolve(nlineq, [1, 1], opt)
nlineq(x)                                            % check
                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          3               5                             6               1
     1          6        0.929784       0.849837           4.99               1
     2          9      0.00256326       0.189768          0.235               1
     3         12     2.83323e-08        0.01091       0.000777               1
     4         15     3.53384e-18    3.64426e-05       8.68e-09               1

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the default value of the function tolerance, and
the problem appears regular as measured by the gradient.

x =

    1.6330    1.1547

ans =

   1.0e-08 *

    0.1329   -0.1329

Data Fitting

Let us see next how to fit a curve into the given data. A good way to find suitable commands is to use a search engine or by typing lookfor fit to list all commands whose description contains the word 'fit'. Out of the listed commands, at least polyfit ja lsqcurvefit seem suitable for curve fitting. Let us generate some data and try these commands out.

x = 5*rand(200, 1);                       % 200 uniformly distributed points in the interval [0,5].
y = 2*x.^2 - 5*x + 1 + 0.5*randn(200, 1); % 2nd order polynomial 2x^2-5x+1 and noise.
pcoef = polyfit(x,y,2)                    % Try fitting 2nd order polynomial obtaining coefficients.
xsorted = sort(x);                        % Order the points for drawing; Matlab draws the connecting
                                          % line segments in the order given by the data points.
ysorted = polyval(pcoef, xsorted);        % Evaluate the polynomial given by the coefficients at the sorted points.
plot(xsorted, ysorted, '-', x, y, '.');
title('Fitting data with polyfit', 'Color', 'k','Fontsize', 15);
plottext=sprintf('Fitted curve: \ny=%3.2fx^{2}%3.2fx+%3.2f', pcoef);
pcoef =

    1.9925   -4.9474    0.9577

Let us do the same fitting with the command lsqcurvefit.

fitf = @(a,x) a(1)*x.^2 + a(2)*x + a(3);
pcoef = lsqcurvefit(fitf, [0 0 0], x, y)
Local minimum found.

Optimization completed because the size of the gradient is less than
the default value of the optimality tolerance.

pcoef =

    1.9925   -4.9474    0.9577

We observe that this command found exactly the same polynomial as polyfit. This is expected as both commands use the method of least squares. Let us next fit a curve a+b*exp(c*x) to new data. The command polyfit cannot be used, but using the command lsqcurvefit is as straightforward as above. Let us also have a look what else the command lsqcurvefit has to offer.

x = 5*rand(200, 1);
y = 2*exp(0.4*x) + 1 + 0.5*randn(200, 1);
fitf = @(a, x) a(1) + a(2)*exp(a(3)*x);
opt = optimset('display', 'iter');
[fcoef,resnorm,res,exitflag,output] = lsqcurvefit(fitf, [0,0,0], x, y, [], [], opt);
xsorted = sort(x);
ysorted = fcoef(1) + fcoef(2)*exp(fcoef(3) * xsorted);
plot(xsorted, ysorted, '-', x, y, '.');
title('Fitting data with lsqcurvefit','Fontsize', 15);
plottext=sprintf('Fitted curve: \ny=%3.2f+%3.2fe^{%3.2fx}', fcoef);
plot(x, res, '.');
title('Residuals','Fontsize', 15);
output.funcCount % the number of evaluations of the function
output.algorithm % the algorithm used
                                         Norm of      First-order 
 Iteration  Func-count     f(x)          step          optimality
     0          4         13443.6                      1.47e+03
     1          8         2616.96        5.20255       3.75e+03      
     2         12         635.566        4.22004       1.06e+03      
     3         16         245.512        1.22828       4.51e+03      
     4         20         83.6733         1.0153       2.53e+03      
     5         24         39.9543       0.085077            104      
     6         28          39.891     0.00629217        0.00452      
     7         32          39.891    3.16258e-05       5.22e-06      

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the default value of the function tolerance.

resnorm =


ans =


exitflag =


ans =


ans =


The command lsqcurvefit can also be used to fit functions of several variables. The convergence could be slow, and you might need to increase the number of iterations.

w = 5*rand(200,1);
z = 5*rand(200,1);
y = 2*exp(0.4*w) + exp(0.2*z) + 1 + 0.1*randn(200,1);
fitf = @(a,x) a(1) + a(2)*exp(a(3)*x(:,1)) + a(4)*exp(a(5)*x(:,2));
fcoef = lsqcurvefit(fitf, [0,0,0,0,0], [w,z], y);
opt = optimset('MaxFunEvals', 30000, 'MaxIter', 5000);
[fcoef,resnorm,res,ef,output] = lsqcurvefit(fitf, [0,0,0,0,0], [w,z], y, [], [], opt);
Local minimum found.

Optimization completed because the size of the gradient is less than
the default value of the optimality tolerance.

Local minimum found.

Optimization completed because the size of the gradient is less than
the default value of the optimality tolerance.

ans =


ans =


When modelling statistical data it may be worthwile to look at functions in Statistics toolbox such as glmfit and mvregress. Let us solve the first example of the section with the function glmfit.

x = 5*rand(200,1);
y = 2*x.^2 - 5*x + 1 + 0.5*randn(200,1);
X = [x, x.^2];                     % Model matrix. Not necessary to write the intercept term.
pcoef = glmfit(X, y, 'normal')
xsorted = sort(x);
ysorted = polyval(flipud(pcoef), xsorted);
plot(xsorted, ysorted, '-', x, y, '.');
title('Fitting data with glmfit', 'Color', 'k','Fontsize', 15);
plottext=sprintf('Fitted curve: \ny=%3.2fx^{2}%3.2fx+%3.2f', flipud(pcoef));
pcoef =


On Interpolation

Interpolation is related to curve fitting. The point of interpolation is to estimate the values of a function between two data points; the interpolated curve will go through all data points unlike in curve fitting. Consider the following data.

x = 1:8; y = exp(0.2*x) + 0.2*transpose(rand(8, 1));

A naive way to find an interpolating function is to fit a polynomial of sufficiently high degree. However, a high-degree polynomial may oscillate significantly between two data points. It is better to use splines.

coef = polyfit(x, y, 7);
xx = 0.6:0.2:8.4;
yy = polyval(coef, xx);
plot(x, y, 'k.', xx, yy, 'b-');
yy = interp1(x, y, xx, 'spline'); % found using lookfor interpolation
plot(xx, yy, 'r-');
legend('data', 'polynomial', 'spline', 'Location', 'NorthWest');
title('Interpolation', 'Color', 'k','Fontsize', 15);
hold off;
Current plot held

Exercise. Do the above spline interpolation using the command fit.