# Lab 1: MATLAB DEVELOPMENT ENVIRONMENT

**Objectives**

Learn about Matlab Development Environment

Learn about the starting and quitting Matlab

Learn about Matlab desktop

Learn about Matlab desktop tools

**Prerequisites**

Students are expected to be familiar with Windows system

Students should have basic knowledge of any text editors

**Requirement**

Matlab software any version

**Remarks**

This lab session will be evaluated and it will carry 2 % of your activity credit points.

All the source codes and executables related to this lab session should be saved in a folder

named Lab1. If not, zero marks will be given. Before you leave the lab make sure you show to your instructor that you have saved your file in Lab1 folder.

**Starting and Quitting MATLAB**

**Starting MATLAB**

On Windows platforms, to start MATLAB, double-click the MATLAB shortcut icon on your Windows desktop.

You can change the directory in which MATLAB starts, define startup options including running a

script upon startup, and reduce startup time in some situations.

**Quitting MATLAB**

To end your MATLAB session, select Exit MATLAB from the File menu in the desktop, or type

quit” in the Command Window. To execute specified functions each time MATLAB quits, such as saving the workspace, you can create and run a name.m” script.

**MATLAB Desktop**

When you start MATLAB, the MATLAB desktop appears, containing tools (graphical user interfaces) for managing files, variables, and applications associated with MATLAB.

You can change the way your desktop looks by opening, closing, moving, and resizing the tools in it.

Use the View menu to open or close the tools. You can also move tools outside the desktop or move

them back into the desktop (docking). All the desktop tools provide common features such as context

menus and keyboard shortcuts.

You can specify certain characteristics for the desktop tools by selecting Preferences from the File

menu. For example, you can specify the font characteristics for Command Window text. For more

information, click the Help button in the Preferences dialog box.

**Desktop Tools**

This section provides an introduction to the MATLAB desktop tools. You can also use MATLAB

functions to perform most of the features found in the desktop tools. The tools are

- Command Window
- Command History
- Start Button and Launch Pad
- Help Browser
- Current Directory Browser
- Workspace Browser
- Array Editor
- Editor/Debugger

**Command Window**

Use the Command Window to enter variables and run functions and M-files. In

this example, MAGIC(4) is called to create an 4-by-4 matrix with the integers from 1 through 4 = 16

with equal row, column, and diagonal sums.

`A=magic(4);`

**Command History**

Statements you enter in the Command Window are logged in the Command History. In the Command

History, you can view previously run statements, and copy and execute selected statements. In green is given the date/time when the execution of the commands took place.

To save the input and output from a MATLAB session to a file, use the diary function as the follows:

`diary filename`

After finishing writing into the file, use

`diary off`

to close the file.

**Start Button**

The MATLAB Start button provides easy access to tools, demos, and documentation. Just click the

button to see the options as shown in Figure below.

**Help Browser**

Use the Help browser to search and view documentation and demos for all your MathWorks products.

The Help browser is a Web browser integrated into the MATLAB desktop that displays HTML

documents.

To open the Help browser, click the help button in the toolbar,

Press F1, or type helpbrowser in

the Command Window and the following image will appear (see Figure below).

The Help browser consists of two panes, the Help Navigator, which you use to find information, and

the display pane, where you view the information. Tabs in the Help Navigator pane provide different

ways to find documentation and demos. Drag the separator bar to adjust the width of the panes. View

documentation in the display pane. Use the close box to hide the pane.

**Help Navigator**

Use the Help Navigator to find information. It includes

- Product filter—set the filter to show documentation only for the products you specify.
- Contents tab—view the titles and tables of contents of documentation for your products.
- Index tab—find specific index entries (selected keywords) in the MathWorks documentation for your products.
- Demos tab—view and run demonstrations for your MathWorks products.
- Search tab—look for a specific word or phrase in the documentation. To get help for a specific function, set the Search type to Function Name.
- Favorites tab—view a list of links to documents you previously designated as favourites.

**Display Pane**

After finding documentation using the Help Navigator, view it in the display pane. While viewing the documentation, you can

- Browse to other pages—Use the arrows at the tops and bottoms of the pages to move through the document, or use the back and forward buttons in the toolbar to go to previously viewed pages.
- Bookmark pages—Click the Add to Favorites button in the toolbar.
- Print pages—Click the print button in the toolbar.
- Find a term in the page—Type a term in the Find in page field in the toolbar and click Go.

Other features available in the display pane are copying information, evaluating a selection, and viewing Web pages.

**For More Help**

In addition to the Help browser, you can use help functions. To get help for a specific function, use doc. For example, doc format displays documentation for the format function in the Help browser.

If you type help followed by the function name, for instance ``help diary”, a briefer form of the documentation appears in the Command Window (see Figure below).

Other means for getting help include contacting Technical Support (http://www.mathworks.com/support) and participating in the newsgroup for MATLAB users, comp.soft-sys.matlab.

**Current Directory Browser**

MATLAB file operations use the current directory and the search path as reference points. Any file you want to run must either be in the current directory or on the search path.

A quick way to view or change the current directory is by using the Current Directory field in the desktop toolbar as shown in figure.

To search for, view, open, and make changes to MATLAB-related directories and files, use the MATLAB Current Directory browser. Alternatively, you can use the functions

`dir`

`cd`

and

`delete`

**Search Path**

MATLAB uses a search path to find M-files and other MATLAB-related files, which are organized in directories on your file system. Any file you want to run in MATLAB must reside in the current directory or in a directory that is on the search path. Add the directories containing files you create to the MATLAB search path. By default, the files supplied with MATLAB and MathWorks toolboxes are included in the search path.

To see which directories are on the search path or to change the search path, select **Search Path** from the **File** menu in the desktop and use the **Set Path** dialog box (shown in the figure below). Alternatively, you can use the path function to view the search path, addpath to add directories to the path, and rmpath to remove the directories from the path. You can have look at the help how to use "path" function in the MATLAB by typing "help path" in the Common Window( see figure below).

**Workspace Browser**

The MATLAB workspace consists of the set of variables (named arrays) built up during a MATLAB session and stored in memory (see figure below). You add variables to the workspace by using functions, running M-files, and loading saved workspaces.

to view workspace and information about each variable, use the Workspace Browser, or use the function "who" and "whos" , for example, typing

`who a*`

in the Command Window will show variable names starting with "a".

To delete variables from the workspace, select the variable and select

**Delete**from the

**Edit**menu. Alternatively, use the clear function.

The workspace is not maintained after you end the MATLAB session. To save the workspace to a file that can be read during a later MATLAB session, select

**Select Workspace As**from the

**File**menu. use the save function. This saves the workspace to a binary file called a MAT-file, which has a .mat extension.

There are options for saving to different formats. To read in a MAT-file, select

**Import Data**from the

**File**menu, or use load function.

**Array Editor**

**Editor/Debugger**

Use the Editor/Debugger to create and debug M-files, which are programs you write to run MATLAB functions. The Editor/Debugger provides a graphical user interface for basic text editing, as well as for M-file debugging (see Figure below).

You can use any text editor to create M-files, such as TextPad, and can use preferences (accessible from the desktop **File** menu) to specify that editor as the default. If you use another editor, you can still use the MATLAB Editor/Debugger for debugging. If you just need to view the contents of an M-file, you can display it in the Command Window by using the type function.

# Lab 2: WORKING WITH MATRICES IN MATLAB

**Objectives:**

Learn about Matrix operations in Matlab

Learn about working with variables, numbers, operators, functions, expressions

Learn about generating matrices, load matrices, create matrices from M-files and concatentation, and delete matrix rows and columns Lear how to use matrices for linear algebra, work with arrays, multivariate data, scalar expansion, and logical subscripting, and use the find function list item

**Prerequisites**

Students are expected to be familiar with Windows system

Students should have basic knowledge of any text editors

**Requirement**

Matlab software any version

**Remarks**

This lab session will be evaluated and it will carry 2 % of your activity credit points.

All the source codes and executables related to this lab session should be saved in a folder named Lab2. If not, zero marks will be given.

Before you leave the lab make sure you show to your instructor that you have saved your file in Lab2 folder.

**Matrices in Matlab**

In MATLAB, a matrix is a two dimensional array of numbers. The 1-by-1 matrices are called scalars, and matrices with only one row or column are called vectors. MATLAB has other ways of storing both numeric and nonnumeric data, but in the beginning, it is usually best to think of everything as a matrix. The operations in MATLAB are designed to be as natural as possible. Where other programming languages work with numbers one at a time, MATLAB allows you to work with entire matrices quickly and easily.

**Entering Matrices**

Start MATLAB and follow along with each example. Start by entering a matrix as a list of its elements. You only have to follow a few basic conventions:

- Separate the elements of a row with blanks or commas.
- Use a semicolon,";", to indicate the end of each row.
- Surround the entire list of elements with square brackets [ ].

To enter square matrix, simply type in the Command Window, for example,

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

And the result should look like (see also figure below):

```
A =
1 3 2
4 2 6
5 9 8
```

Once you have entered the matrix, it is automatically remembered in the MATLAB workspace. You can refer to it simply as A. Now that you have A in the workspace, take a look at what makes it so interesting.

**sum, transpose, and diag**

The first statement to try is

`>> sum(A)`

MATLAB replies with

```
ans =
10 14 16
```

When you do not specify an output variable, MATLAB uses the variable ans, short for answer, to store the results of a calculation. You have computed a row vector containing the sums of the columns of A.

How about the row sums? MATLAB has a preference for working with the columns of a matrix, so the easiest way to get the row sums is to transpose first the matrix, compute the column sums of the transpose, and then transpose the result. The transpose operation is denoted by an apostrophe or single quote as:

`>> A’`

It flips a matrix about its main diagonal and it turns a row vector into a column vector:

```
ans =
1 4 5
3 2 9
2 6 8
```

Then,

`>> sum(A’)’`

will produce:

```
ans =
6
12
22
```

You can check it easily, that this is the sum along the rows of the matrix A.

The sum of the elements on the main diagonal is obtained with the sum and the diag functions.

`>> diag(A)`

produces

```
ans =
1
2
8
```

and

`>> sum(diag(A))`

produces

```
ans =
11
```

**Subscripts**

The element in row i and column j of A is denoted by A(i,j). For example, A(1,2) is the number in the first row and second column. For our square matrix, A(1,2) is 3, which can be verified by typing:

`>> A(1,2)`

```
ans =
3
```

So to compute the sum of the elements in the first column of A, type:

`>> A(1,1) + A(2,1) + A(3,1)`

which produces:

```
ans =
10
```

but is not the most elegant way of summing a single column.

It is also possible to refer to the elements of a matrix with a single subscript, A(k). This is the usual way of referencing row and column vectors. But it can also apply to a fully two-dimensional matrix, in which case the array is regarded as one long column vector formed from the columns of the original matrix. So, for our matrix, A(4) is another way of referring to the value

3 stored in A(1,2):

```
>> A(1,2)
ans =
3
```

```
>> A(4)
ans =
3
```

If you try to use the value of an element outside of the matrix, it is an error.

`>> A(4,5)`

Index exceeds matrix dimensions.

**The Colon Operator**

The colon, :, is one of the most important MATLAB operators. It occurs in several different forms. The expression

1:5 is a row vector containing the integers from 1 to 5

```
ans =
1 2 3 4
```

5

To obtain non unit spacing, specify an increment. For example,

```
>> 10:5:30
ans =
10 15 20 25
```

30

and

`>>0:pi/4:pi`

is

`0 0.7854 1.5708 2.3562 3.1416`

Subscript expressions involving colons refer to portions of a matrix. A(1:k,j) is the first k elements of the jth column of A. So

`>>sum(A(1:3,3))`

Computes the sum of the third column:

```
ans =
16
```

But there is a better way. The colon by itself refers to all the elements in a row or column of a matrix and the keyword end refers to the last row or column. So

`>>sum(A(:,end))`

computes the sum of the elements in the last column of A.

```
ans =
16
```

**Expressions**

Like most other programming languages, MATLAB provides mathematical expressions, but unlike most programming languages, these expressions involve entire matrices. The building blocks of expressions are

- Variables
- Numbers
- Operators
- Functions

**Variables**

MATLAB does not require any type declarations or dimension statements. When MATLAB encounters a new variable name, it automatically creates the variable and allocates the appropriate amount of storage. If the variable already exists, MATLAB changes its contents and, if necessary, allocates new storage. For example,

`>> i = 25`

creates a 1-by-1 matrix named i and stores the value 25 in its single element.

Variable names consist of a letter, followed by any number of letters, digits, or underscores. MATLAB uses only the first 31 characters of a variable name. MATLAB is case sensitive; it distinguishes between uppercase and lowercase letters. A and a are not the same variable. To view the matrix assigned to any variable, simply enter the variable name.

**Numbers**

MATLAB uses conventional decimal notation, with an optional decimal point and leading plus or minus sign, for numbers. Scientific notation uses the letter e to specify a power-of-ten scale factor. Imaginary numbers use either i or j as a suffix.

Some examples of legal numbers are:

3 -99 0.0001 9.6397238 1.60210e-20 6.02252e23 1i -3.14159j 3e5i

**Operators**

Expressions use familiar arithmetic operators and precedence rules.

(+) Addition

(-) Subtraction

(*) Multiplication

(/) Division

(\) Left division (described in “Matrices and Linear Algebra” in the MATLAB documentation)

(^) Power

(')Complex conjugate transpose

( ) Specify evaluation order

**Functions**

MATLAB provides a large number of standard elementary mathematical functions, including abs, sqrt, exp, and sin. Taking the square root or logarithm of a negative number is not an error; the appropriate complex result is produced automatically. MATLAB also provides many more advanced mathematical functions, including Bessel and gamma functions. Most of these functions accept complex arguments. For a list of the elementary mathematical functions, type

help elfun

For a list of more advanced mathematical and matrix functions, type

```
help specfun
help elmat
```

Some of the functions, like sqrt and sin, are *built in*. They are part of the MATLAB core so they are very efficient, but the computational details are not readily accessible. Other functions, like gamma and sinh, are implemented in M-files. You can see the code and even modify it if you want.

Several special functions provide values of useful constants.

pi 3.14159265…

i Imaginary unit, √-1

j Same as i

eps Floating-point relative precision, 2-52

realmin Smallest floating-point number, 2-1022

realmax Largest floating-point number, (2-ε)21023

Inf Infinity

NaN Not-a-number

Infinity is generated by dividing a nonzero value by zero, or by evaluating well defined mathematical expressions that *overflow,* i.e., exceed realmax. Not-a-number is generated by trying to evaluate expressions like 0/0 or Inf-Inf that do not have well defined mathematical values.

The function names are not reserved. It is possible to overwrite any of them with a new variable, such as

`>>eps = 1.e-6`

and then use that value in subsequent calculations. The original function can be restored with

`>>clear eps`

**Examples of Expressions**

You have already seen several examples of MATLAB expressions. Here are a

few more examples, and the resulting values.

```
>> r=1+sqrt(2)/5
r =
1.2828
```

**Working with matrices**

**Generating Matrices**

zeros All zeros

ones All ones

rand Uniformly distributed random elements

randn Normally distributed random elements

```
>> z=zeros(3,3)
z =
0 0 0
0 0 0
0 0 0
```

```
>> T=2*ones(2,3)
T =
2 2 2
2 2 2
```

```
>> N=rand(1,5)
N =
0.8147 0.9058 0.1270 0.9134 0.6324
```

```
>> G=randn(1,5)
G =
-1.3077 -0.4336 0.3426 3.5784 2.7694
```

**M-Files**

You can create your own matrices using M-files, which are text files containing MATLAB code. Use the MATLAB Editor or another text editor to create a file name it matFile.m.

For example, create a file containing these five lines.

```
A = [ ...
16.0 3.0 2.0 13.0
5.0 10.0 11.0 8.0
9.0 6.0 7.0 12.0
4.0 15.0 14.0 1.0 ];
```

Then the statement

matFile

reads the file and creates a variable, A, containing our example matrix.

```
A =
16 3 2 13
5 10 11 8
9 6 7 12
4 15 14 1
```

**Linear Algebra**

The mathematical operations defined on matrices are the subject of linear algebra. Adding a matrix to its transpose produces a symmetric matrix.

```
A =
16 3 2 13
5 10 11 8
9 6 7 12
4 15 14 1
```

```
>> A+A'
ans =
32 8 11 17
8 20 17 23
11 17 14 26
17 23 26 2
```

The multiplication symbol, *, denotes the matrix multiplication involving inner products between rows and columns. Multiplying the transpose of a matrix by the original matrix also produces a symmetric matrix.

```
>> A'*A
ans =
378 212 206 360
212 370 368 206
206 368 370 212
360 206 212 378
```

The determinant of this particular matrix happens to be zero, indicating that

the matrix is singular.

```
>> det(A)
ans =
1.0871e-012
```

Since the matrix is singular, it does not have an inverse. If you try to compute the inverse with you will get a warning message

`>> inv(A)`

Warning: Matrix is close to singular or badly scaled.

Results may be inaccurate. RCOND = 9.796086e-018.

```
ans =
1.0e+015 *
0.1251 0.3753 -0.3753 -0.1251
-0.3753 -1.1259 1.1259 0.3753
0.3753 1.1259 -1.1259 -0.3753
-0.1251 -0.3753 0.3753 0.1251
```

The eigenvalues of the magic square are interesting.

```
>> eig(A)
ans =
34.0000
8.0000
0.0000
-8.0000
```

If the square matrix is multiplied by itself with array multiplication A.*A the result is an array containing the squares of the integers from 1 to 16, in an unusual order.

```
>> A.*A
ans =
256 9 4 169
25 100 121 64
81 36 49 144
16 225 196 1
```

# Lab 3: WORKING WITH MATRICES IN MATLAB

**Objectives:**

Learn more about Matrix operations in Matlab

Learn about creating plots in Matlab.

**Prerequisites**

Students are expected to be familiar with Windows system

Students should have basic knowledge of any text editors

**Requirement**

Matlab software any version

**Remarks**

This lab session will be evaluated and it will carry 2 % of your activity credit points.

All the source codes and executables related to this lab session should be saved in a folder named Lab3. If not, zero marks will be given.

Before you leave the lab make sure you show to your instructor that you have saved your file in Lab3 folder.

**More about matrices in Matlab**

We can reshape matrices in Matlab, for example we will consider reshaping a 3-by-4 matrix A to have dimensions 2-by-6:

`>> A = rand(3,4);`

```
>> A
A =
0.8944 0.9274 0.6183 0.1248
0.1375 0.9175 0.3433 0.7306
0.3900 0.7136 0.9360 0.6465
```

```
>> B=reshape(A, 2, 6)
B =
0.8944 0.3900 0.9175 0.6183 0.9360 0.7306
0.1375 0.9274 0.7136 0.3433 0.1248 0.6465
```

Transpose A so that the row elements become columns. You can use either the transpose function or the transpose operator (.') to do this:

```
>> A'
ans =
0.8944 0.1375 0.3900
0.9274 0.9175 0.7136
0.6183 0.3433 0.9360
0.1248 0.7306 0.6465
```

Rotate the matrix by 90 degrees:

```
>> B=rot90(A)
B =
0.1248 0.7306 0.6465
0.6183 0.3433 0.9360
0.9274 0.9175 0.7136
0.8944 0.1375 0.3900
```

Flip A in a left-to-right direction:

```
>> B=fliplr(A)
B =
0.1248 0.6183 0.9274 0.8944
0.7306 0.3433 0.9175 0.1375
0.6465 0.9360 0.7136 0.3900
```

TF = isfinite(A) returns an array the same size as A containing logical 1 (true) where the elements of the array A are finite and logical 0 (false) where they are infinite or NaN.

```
>> TF = isfinite(A)
TF =
1 1 1 1
1 1 1 1
1 1 1 1
```

d = size(X) returns the sizes of each dimension of array A in a vector d with ndims(A) elements. If A is a scalar, which MATLAB software regards as a 1-by-1 array, size(A) returns the vector [1 1]. [m,n] = size(A) returns the size of matrix A in separate variables m and n:

```
>> [m,n]=size(A)
m =
3
n =
4
```

[V,D] = eig(A) produces matrices of eigenvalues (D) and eigenvectors (V) of a square matrix A, so that A*V = V*D. For example:

```
>> A=rand(3,3)
A =
0.8332 0.8352 0.9791
0.3983 0.3225 0.5493
0.7498 0.5523 0.3304
```

```
>> [V,D]=eig(A)
V =
-0.7744 -0.4620 -0.3334
-0.3768 0.8400 -0.4880
-0.5083 -0.2845 0.8067
D =
1.8822 0 0
0 -0.0826 0
0 0 -0.3135
```

[row,col,v] = find(X) returns a column or row vector v of the nonzero entries in X, as well as row and column indices. If X is a logical expression, then v is a logical array. Output v contains the non-zero elements of the logical array obtained by evaluating the expression X. For example:

```
>> A=magic(3)
A =
8 1 6
3 5 7
4 9 2
```

```
>> [r,c,v]= find(A>5);
>> r', c', v'
ans =
1 3 1 2
ans =
1 2 3 3
ans =
1 1 1 1
```

The scalar product of the vectors A and B is computed as:

dot(A,B)

where A and B must be vectors of the same length. When A and B are both column vectors, DOT(A,B) is the same as A'*B:

```
>> A=rand(5,1)
A =
0.6947
0.9727
0.3278
0.8378
0.7391
```

```
>> B=rand(5,1)
B =
0.9542
0.0319
0.3569
0.6627
0.2815
```

```
>> C=dot(A,B)
C =
1.5741
```

```
>> D=A'*B
D =
1.5741
```

If A and B are row vectors, then dot(A,B) is equivalent as: A*B’. For example,

```
>> A=rand(1,5)
A =
0.2304 0.7111 0.6246 0.5906 0.6604
```

```
>> B=rand(1,5)
B =
0.0476 0.3488 0.4513 0.2409 0.7150
```

```
>> C=dot(A,B)
C =
1.1554
```

```
>> D=A*B'
D =
1.1554
```

We can also find the inverse of a matrix. You must be careful, however, since the operations are numerical manipulations done on digital computers. In the example, the matrix A is not a full matrix, but matlab's inverse routine will still return a matrix.

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

```
>> inv(A)
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 3.469447e-018.
ans =
1.0e+015 *
-2.7022 4.5036 -1.8014
5.4043 -9.0072 3.6029
-2.7022 4.5036 -1.8014
```

By the way, Matlab is case sensitive. This is another potential source of problems when you start building complicated algorithms.

```
>> inv(a)
??? Undefined function or variable a.
```

There are also routines that let you find solutions to equations. For example, if Ax=b and you want to find x, a slow way to find x is to simply invert A and perform a left multiply on both sides (more on that later). It turns out that there are more efficient and more stable methods to do this (L/U decomposition with pivoting, for example). Matlab has special commands that will do this for you.

Before finding the approximations to linear systems, it is important to remember that if A and B are both matrices, then AB is not necessarily equal to BA. To distinguish the difference between solving systems that have a right or left multiply, Matlab uses two different operators, "/" and "\". Examples of their use are given below. It is left as an exercise for you to figure out which one is doing what.

```
>> v = [1 3 5]'
v =
1
3
5
```

```
>> x = A\v
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 4.565062e-18
```

```
x =
1.0e+15 *
1.8014
-3.6029
1.8014
```

```
>> x = B\v
x =
2
1
-1
```

```
>> B*x
ans =
1
3
5
```

```
>> x1 = v'/B
x1 =
4.0000 -3.0000 1.0000
>> x1*B
ans =
1.0000 3.0000 5.0000
```

**Vectors and Matrices and Matlab Functions**

If you pass a vector to a predefined math function, it will return a vector of the same size, and each entry is found by performing the specified operation on the corresponding entry of the original vector: (Notice that the second command has a ";" at the end of the line. This tells Matlab that it should not print out the result.)

```
>> v = [1 2 3]';
>> sin(v)
ans =
0.8415
0.9093
0.1411
```

```
>> log(v)
ans =
0
0.6931
1.0986
```

The ability to work with these vector functions is one of the advantages of Matlab. Now complex operations can be defined that can be done quickly and easily. In the following example a very large vector is defined and can be easily manipulated.

```
>> x = [0:0.1:100]
x =
Columns 1 through 7
0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000
[stuff deleted]
Columns 995 through 1001
99.4000 99.5000 99.6000 99.7000 99.8000 99.9000 100.0000
```

`>> y = sin(x).*x./(1+cos(x));`

**Creating a Plot**

The plot function has different forms, depending on the input arguments. If y is a vector, plot(y) produces a piecewise linear graph of the elements of y versus the index of the elements of y. If you specify two vectors as arguments, plot(x,y) produces a graph of y versus x.

For example, these statements use the colon operator to create a vector of x values ranging from zero to compute the cosine of these values, and plot the result.

```
>> x=0:0.1:(2*pi);
>> y=cos(x);
>> plot(x,y)
```

We can also label the axes and add a title as the following:

```
>> xlabel('values of x');
>> ylabel('cos(x)');
>> title('cosine function plot')
```

You can also specify the font size as follows:

```
>> xlabel('values of x','FontSize',12);
>> ylabel('cos(x)', 'FontSize',12);
>> title('cosine function plot','FontSize',12)
```

**Multiple Data Sets in One Graph**

Multiple x-y pair arguments create multiple graphs with a single call to plot. MATLAB automatically cycles through a predefined (but user settable) list of colors to allow discrimination among sets of data. For example, these statements plot three related functions of x, each curve in a separate distinguishing color.

```
>> y2 = cos(5*x);
>> y3 = cos(2.5*x);
>> plot(x,y,x,y2,x,y3)
```

The legend command provides an easy way to identify the individual plots (see also Figure 2).

`>> legend('cos(x)','cos(5*x)','cos(2.5*x)')`

**Specifying Line Styles and Colors**

It is possible to specify color, line styles, and markers (such as plus signs or circles) when you plot your data using the plot command.

`plot(x,y,'color_style_marker')`

color_style_marker is a string containing from one to four characters (enclosed in single quotation marks) constructed from a color, a line style, and a marker type:

Color strings are 'c', 'm', 'y', 'r', 'g', 'b', 'w', and 'k'. These correspond to cyan, magenta, yellow, red, green, blue, white, and black.

Linestyle strings are '-' for solid, '—' for dashed, ':' for dotted, '-.' For dash-dot. Omit the linestyle for no line.

The marker types are '+', 'o', '*', and 'x' and the filled marker types are's' for square, 'd' for diamond, '^' for up triangle, 'v' for down triangle, '>' for right triangle, '<' for left triangle, 'p' for pentagram, 'h' for hexagram, and none for no marker.

**Plotting Lines and Markers**

If you specify a marker type but not a linestyle, MATLAB draws only the marker. For example,

`plot(x,y,'ko')`

plots black circle at each data point, but does not connect the markers with a line. The statement

`plot(x,y,'b-o')`

plots a blue solid line and places circle sign markers at each data point. You may want to use fewer data points to plot the markers than you use to plot the lines. This example plots the data twice using a different number of points for the solid line and marker plots (see Figure below).

```
x1 = 0:pi/100:2*pi;
y1=cos(x1);
x2 = 0:pi/10:2*pi;
y2=cos(x2);
plot(x1,y1,'r-',x2,y2,'ro')
```

**Adding Plots to an Existing Graph**

The hold command enables you to add more than one plot to an existing graph. When you type

`hold on`

after each plot command MATLAB does not replace the existing graph, but it adds the new data to the current graph, rescaling the axes if needed.

For example the following Matlab commands will create the same plot as in previous figure.

```
x1 = 0:pi/100:2*pi;
y1=cos(x1);
plot(x1, y1, ‘r-‘);
hold on;
x2 = 0:pi/10:2*pi;
y2=cos(x2);
plot(x2,y2,'ro')
hold on;
```

# Lab4: WORKING WITH PLOTS AND CODING IN MATLAB

**Objectives:**

Learn more about creating plots in Matlab.

Learn about subplot.

Learn about coding in Matlab.

**Prerequisites**

Students are expected to be familiar with Windows system

Students should have basic knowledge of any text editors

Requirement

**Matlab software any version**

Remarks

This lab session will be evaluated and it will carry 2 % of your activity credit points.

All the source codes and executables related to this lab session should be saved in a folder named Lab4. If not, zero marks will be given.

Before you leave the lab make sure you show to your instructor that you have saved your file in Lab4 folder.

**More about plots in Matlab**

**Figure Windows**

If there are no figure windows already on the screen, plotting functions automatically open a new figure window. If a figure window exists, MATLAB uses that window for graphics output. If there are multiple figure windows open, MATLAB targets the one that is designated the “current figure” (the last figure used or clicked in).

To make an existing figure window the current figure, you can click the mouse while the pointer is in that window or you can type

`figure(n)`

where n is the number in the figure title bar. The results of subsequent graphics commands are displayed in this window.

To open a new figure window and make it the current figure, type

`figure`

Multiple Plots in One Figure

The subplot command enables one to display multiple plots in the same figure or print them on the same piece of paper. Typing

subplot(m,n,p)

partitions the figure window into an m-by-n matrix of small subplots and selects the pth subplot for the current plot. The axes are counted along the top row of the figure window, then the second row, etc. For example,

```
income = [3.2 4.1 5.0 5.6];
outgo = [2.5 4.0 3.35 4.9];
subplot(2,1,1); plot(income)
title('Income')
subplot(2,1,2); plot(outgo)
title('Outgo')
```

Another example is given below creating 4 subplots (see also Figure 2):

```
t= 0:pi/10:2*pi;
[X,Y,Z] = cylinder(4*cos(t));
subplot(2,2,1); mesh(X)
subplot(2,2,2); mesh(Y)
subplot(2,2,3); mesh(Z)
subplot(2,2,4); mesh(X,Y,Z)
```

**Controlling the Axes**

The axis command supports a number of options for setting the scaling, orientation, and aspect ratio of plots.

**Setting Axis Limits**

The axis command enables you to specify the axes limits, for example:

`axis([xmin xmax ymin ymax])`

sets the limits for the x- and y-axis of the current axes. For three-dimensional plot, we may use:

`axis([xmin xmax ymin ymax zmin zmax])`

While the command

`axis auto`

enables MATLAB automatic limit selection.

**Setting Axis Aspect Ratio**

axis also enables you to specify a number of predefined modes. For example,

`axis square`

makes the x-axes and y-axes the same length.

`axis equal`

makes the individual tick mark increments on the x- and y-axes the same length. This means

`plot(exp(i*[0:pi/10:2*pi]))`

followed by either axis square or axis equal turns the oval into a proper circle.

`axis auto normal`

returns the axis scaling to its default, automatic mode.

**Setting Axis Visibility**

You can use the axis command to make the axis visible or invisible.

`axis on`

makes the axis visible. This is the default.

`axis off`

makes the axis invisible.

**Setting Grid Lines**

The grid command toggles grid lines on and off. The statement

`grid on`

turns the grid lines on and

`grid off`

turns them back off again.

**Axis Labels and Titles**

The xlabel, ylabel, and zlabel commands add x-, y-, and z-axis labels, respectively. The title command adds a title at the top of the figure and the text function inserts text anywhere in the figure. A subset of TeX notation produces Greek letters. You can also set these options interactively.

For example,

```
t = -pi:pi/100:pi;
y = sin(t).*cos(t);
plot(t,y)
axis([-pi pi -1 1])
xlabel('-\pi \leq {\itt} \leq \pi')
ylabel('sin(t)*cos(t)')
title('Graph of the sine \times cosine functions')
text(-2.5,0.6,'{\itMax}')
text(-0.9,-0.6,'{\itMin}')
text(0.9,0.6,'{\itMax}')
text(2.5,-0.6,'{\itMin}')
```

You can save the figure in different formats. To do this, select **Save** from the **File** menu. To save it using a graphics format, such as JPEG, PNG, etc., for use with other applications, select **Export** from the **File** menu, then graph shown in figure will appear. Select **Export** to save the file using the desired format.

**Programming with Matlab**

**Flow Control**

MATLAB has several flow control constructs:

• “if”

• “switch and case”

• “for”

• while

• “continue”

• “break”

**if**

The if statement evaluates a logical expression and executes a group of statements when the expression is true. The optional elseif and else keywords provide for the execution of alternate groups of statements. An end keyword, which matches the if, terminates the last group of statements. The groups of statements are delineated by the four keywords—no braces or brackets are involved.

The MATLAB algorithm for generating a random square n-by-n involves three different cases: when n is odd, when n is even but not divisible by 4, or when n is divisible by 4. This is described by

```
if rem(n,2) ~= 0
M = randn(n,n);
elseif rem(n,4) ~= 0
M = ones(n,n)
else
M = zeros(n,n)
End
```

In this example, the three cases are mutually exclusive, but if they weren’t, the first true condition would be executed. REM function is the remainder after division.

**switch and case**

The switch statement executes groups of statements based on the value of a variable or expression. The keywords case and otherwise delineate the groups. Only the first matching case is executed. There must always be an end to match the switch.

The previous example can also be written as:

```
switch (rem(n,4)==0) + (rem(n,2)==0)
case 0
M = randn(n,n)
case 1
M = ones(n,n)
case 2
M = zeros(n,n)
otherwise
error('This is impossible')
end
```

**for**

The for loop repeats a group of statements a fixed, predetermined number of times. A matching end delineates the statements.

```
for n = 3:32
r(n) = randn(1,1);
end
r
```

The semicolon terminating the inner statement suppresses repeated printing, and the r after the loop displays the final result. It is a good idea to indent the loops for readability, especially when they are nested. For this, use Text →Smart Indent.

```
for i = 1:m
for j = 1:n
H(i,j) = 1/(i+j);
end
end
```

**while**

The while loop repeats a group of statements an indefinite number of times under control of a logical condition. A matching end delineates the statements.

```
i=0;
n=20;
s=0;
while i < n
i=i+1;
s = s+i^2;
end
s
```

**continue**

The continue statement passes control to the next iteration of the for or while loop in which it appears, skipping any remaining statements in the body of the loop. In nested loops, continue passes control to the next iteration of the for or while loop enclosing it.

The example below shows a continue loop that counts the lines of code in the file, magic.m, skipping all blank lines and comments. A continue statement is used to advance to the next line in magic.m without incrementing the count whenever a blank line or comment line is encountered.

```
fid = fopen('splots6.m','r');
count = 0;
while ~feof(fid)
line = fgetl(fid);
if isempty(line) | strncmp(line,'%',1)
continue
end
count = count + 1;
end
disp(sprintf('%d lines',count));
```

where ``splots6.m” looks as follows:

```
i=0;
n=20;
s=0;
while i < n
i=i+1;
s = s+i^2;
end
s
```

**break**

The break statement lets you exit early from a for or while loop. In nested loops, break exits from the innermost loop only. Here is an improvement on the example from the previous section. Why is this use of break a good idea?

```
i=0;
n=20;
s=0;
while i < n
i=i+1;
s = s+i^2;
if s>= 100
break;
end
end
s
```

# Lab 5: SYSTEM MODELING AND NUMERICAL INTEGRATION IN MATLAB

**Objectives:**

Learn more about modelling a simple system in matlab.

Learn about Euler method of numerical integration in matlab.

**Prerequisites**

Students are expected to be familiar with Windows system

Students should have basic knowledge of any text editors

**Requirement**

Matlab software any version

**Remarks**

This lab session will be evaluated and it will carry 2 % of your activity credit points.

All the source codes and executables related to this lab session should be saved in a folder named Lab5. If not, zero marks will be given.

Before you leave the lab make sure you show to your instructor that you have saved your file in Lab5 folder.

**Example 1: Conveyor System**

Consider a factory conveyor system in which boxes arrive at the rate of one box each Δt. Each box has a weight of w_{1}, w_{2} , w_{3} . However, there are twice as many w_{1} boxes and w_{3} boxes as w_{2} boxes. How do we model and simulate this?

The Probability mass distribution is:

(1)where W is a weight random variable that can take on one of the three discrete values. The notation P(W=w) gives the probability that the random variable W is w.

The set {w_{1} w_{2} w_{3} } is called the sample space of , and is the set of all possible weights. According to the description, these boxes arrive every Δt seconds, so t = Δt k gives the continuous time measured in successive k-values, assuming the initial time is zero. However, how do we describe the system output?

An interesting problem would be the weight W(t) of the boxes as they arrive. W(t) is a non-deterministic variable, and we can only hope to simulate the behavior of this variable. If we denote with W(k) the weight of the k-th event. This can be accomplished by using the random (RAND) function, which is a random number generator that provides uniformly random distributed variates such that

(2)```
function lab5_ex1
clear all;
N=100;
w=zeros(1,N);
w1=1; w2=2; w3=3;
Nw=zeros(1,3);
P=zeros(1,3);
for k=1:N
r=10*rand(1,1);
if r < 4
w(k) = w1;
Nw(1)=Nw(1)+1;
elseif r >= 4 && r < 6
w(k) = w2;
Nw(2)=Nw(2)+1;
else
w(k) = w3;
Nw(3)=Nw(3)+1;
end
end
P=Nw/N;
Pw=[0.4 0.2 0.4];
figure(1);
hold on;
subplot(2,1,1); bar(P,'grouped','r');
subplot(2,1,2); bar(Pw,'grouped','b');
colormap jet;
figure(2);
hold on;
stairs(cumsum(w), 'k-');
```

**Euler Numerical integration method:**

A matlab code for numerical integration using Euler 1st order method is shown below:

```
function [wi, ti] = euler ( RHS, t0, x0, tf, N )
%EULER approximate the solution of the initial value problem
%
% x'(t) = RHS( t, x ), x(t0) = x0
%
% using Euler's method - this routine will work for a system
% of first-order equations as well as for a single equation
%
% calling sequences:
% [wi, ti] = euler ( RHS, t0, x0, tf, N )
% euler ( RHS, t0, x0, tf, N )
%
% inputs:
% RHS string containing name of m-file defining the
% right-hand side of the differential equation; the
% m-file must take two inputs - first, the value of
% the independent variable; second, the value of the
% dependent variable
% t0 initial value of the independent variable
% x0 initial value of the dependent variable(s)
% if solving a system of equations, this should be a
% row vector containing all initial values
% tf final value of the independent variable
% N number of uniformly sized time steps to be taken to
% advance the solution from t = t0 to t = tf
%
% output:
% wi vector / matrix containing values of the approximate
% solution to the differential equation
% ti vector containing the values of the independent
% variable at which an approximate solution has been
% obtained
%
neqn = length ( x0 );
ti = linspace ( t0, tf, N+1 );
wi = [ zeros( neqn, N+1 ) ];
wi(1:neqn, 1) = x0';
h = ( tf - t0 ) / N;
for i = 1:N
x0 = x0 + h * feval ( RHS, t0, x0 );
t0 = t0 + h;
wi(1:neqn,i+1) = x0';
end;
```

**Example**

Consider the following differential equation:

(3)With initial conditions:

(4)A matlab function for returning a value of $\frac{dx}{dt}$ is given in the following:

```
function f = example(t,x)
f=t*(x^2);
end
```

A matlab testing Euler method for different values of the time step is given below:

```
function Eulerdemo
x0=3;
t0=1;
tf=1.3;
h=0.02;
N=floor((tf-t0)/h)
[x4, t4]=euler('example', t0, x0, tf, N);
xt = 6./(5-3*t4.^2);
plot(t4,x4, 'r-x', 'LineWidth',2,...
'MarkerEdgeColor','r',...
'MarkerFaceColor','r',...
'MarkerSize',10);
hold on;
plot(t4, xt, 'k-', 'LineWidth',2);
hold on;
h=0.04;
N=floor((tf-t0)/h)
[x4, t4]=euler('example', t0, x0, tf, N);
plot(t4,x4, 'b-o', 'LineWidth',2,...
'MarkerEdgeColor','r',...
'MarkerFaceColor','r',...
'MarkerSize',10);
hold on;
h=0.08;
N=floor((tf-t0)/h)
[x4, t4]=euler('example', t0, x0, tf, N);
plot(t4,x4, 'g-s', 'LineWidth',2,...
'MarkerEdgeColor','r',...
'MarkerFaceColor','r',...
'MarkerSize',10);
hold on;
xlim([1 1.25]);
```

Graphically results are presented in figure:

The analytical solution of the above differential equation is:

(5)**+ Lab 6: TAYLOR AND RUNGE-KUTTA NUMERICAL INTEGRATORS IN MATLAB**

**Objectives:**

- Learn more about numerical integration in matlab.
- Comparison of numerical integrators in matlab.

**Prerequisites**

- Students are expected to be familiar with Windows system
- Students should have basic knowledge of any text editors

**Requirement**

Matlab software any version

**Remarks**

This lab session will be evaluated and it will carry 2 % of your activity credit points.

All the source codes and executables related to this lab session should be saved in a

folder named Lab6. If not, zero marks will be given.

Before you leave the lab make sure you show to your instructor that you have saved

your file in Lab6 folder.

**Example 1:**

Consider again the following differential equation:

(6)with initial condition:

(7)A matlab function for returning a value of $\frac{dx}{dt}$ is given in the following:

```
function f = example(t,x)
f=t*(x^2);
end
```

**Euler Numerical integration method:**

A matlab code for numerical integration using Euler 1st order method is shown below:

```
function [wi, ti] = euler ( RHS, t0, x0, tf, N )
%EULER approximate the solution of the initial value problem
%
% x'(t) = RHS( t, x ), x(t0) = x0
%
% using Euler's method - this routine will work for a system
% of first-order equations as well as for a single equation
%
% calling sequences:
% [wi, ti] = euler ( RHS, t0, x0, tf, N )
% euler ( RHS, t0, x0, tf, N )
%
% inputs:
% RHS string containing name of m-file defining the
% right-hand side of the differential equation; the
% m-file must take two inputs - first, the value of
% the independent variable; second, the value of the
% dependent variable
% t0 initial value of the independent variable
% x0 initial value of the dependent variable(s)
% if solving a system of equations, this should be a
% row vector containing all initial values
% tf final value of the independent variable
% N number of uniformly sized time steps to be taken to
% advance the solution from t = t0 to t = tf
%
% output:
% wi vector / matrix containing values of the approximate
% solution to the differential equation
% ti vector containing the values of the independent
% variable at which an approximate solution has been
% obtained
%
neqn = length ( x0 );
ti = linspace ( t0, tf, N+1 );
wi = [ zeros( neqn, N+1 ) ];
wi(1:neqn, 1) = x0';
h = ( tf - t0 ) / N;
for i = 1:N
x0 = x0 + h * feval ( RHS, t0, x0 );
t0 = t0 + h;
wi(1:neqn,i+1) = x0';
end;
```

```
Taylor 2nd order numerical method in matlab:
function [wi, ti] = taylor2nd ( RHS, t0, x0, tf, N )
%TAYLOR2ND approximate the solution of the initial value problem
% S
% x'(t) = RHS( t, x ), x(t0) = x0
%
% using the second-order Taylor method - this routine will
% work for a system of first-order equations as well as for
% a single equation
%
% calling sequences:
% [wi, ti] = taylor2nd (RHS, t0, x0, tf, N )
% taylor2nd ( RHS, t0, x0, tf, N )
%
% inputs:
% RHS string containing name of m-file defining the
% right-hand side of the differential equation and
% its derivative with respect to the independent
% variable;
% the prototype for the m-file should be
% [f, ft] = RHS ( t, x )
% where f is the value of the right-hand side function
% and ft is the value of the derivative with respect
% to the independent variable
% t0 initial value of the independent variable
% x0 initial value of the dependent variable(s)
% if solving a system of equations, this should be a
% row vector containing all initial values
% tf final value of the independent variable
% N number of uniformly sized time steps to be taken to
% advance the solution from t = t0 to t = tf
%
% output:
% wi vector / matrix containing values of the approximate
solution to the differential equation
% ti vector containing the values of the independent
% variable at which an approximate solution has been
% obtained
%
neqn = length ( x0 );
ti = linspace ( t0, tf, N+1 );
wi = [ zeros( neqn, N+1 ) ];
wi(1:neqn, 1) = x0';
h = ( tf - t0 ) / N;
for i = 1:N
[f ft] = feval ( RHS, t0, x0 );
x0 = x0 + h * f + h^2 * ft / 2;
t0 = t0 + h;
wi(1:neqn,i+1) = x0';
end;
```

A matlab testing Euler method for different values of the time step is given below:

```
function Eulerdemo
clear all;
x0=3;
t0=1;
tf=1.3;
h=0.02;
N=floor((tf-t0)/h)
[x4, t4]=euler('example', t0, x0, tf, N);
xt = 6./(5-3*t4.^2);
plot(t4,x4, 'r-x', 'LineWidth',2,...
'MarkerEdgeColor','r',...
'MarkerFaceColor','r',...
'MarkerSize',10);
hold on;
plot(t4, xt, 'k-', 'LineWidth',2);
hold on;
h=0.04;
N=floor((tf-t0)/h)
[x4, t4]=euler('example', t0, x0, tf, N);
plot(t4,x4, 'b-o', 'LineWidth',2,...
'MarkerEdgeColor','r',...
'MarkerFaceColor','r',...
'MarkerSize',10);
hold on;
h=0.08;
N=floor((tf-t0)/h)
[x4, t4]=euler('example', t0, x0, tf, N);
hold on;
xlim([1 1.25]);
```

**Runge-Kutta 4th order:**

```
function [wi, ti] = rk4 ( RHS, t0, x0, tf, N )
%RK4 approximate the solution of the initial value problem
%
% x'(t) = RHS( t, x ), x(t0) = x0
%
% using the classical fourth-order Runge-Kutta method - this
% routine will work for a system of first-order equations as
% well as for a single equation
%
% calling sequences:
plot(t4,x4, 'g-s', 'LineWidth',2,...
'MarkerEdgeColor','r',...
'MarkerFaceColor','r',...
'MarkerSize',10);
% [wi, ti] = rk4 ( RHS, t0, x0, tf, N )
% rk4 ( RHS, t0, x0, tf, N )
%
% inputs:
% RHS string containing name of m-file defining the
% right-hand side of the differential equation;
the
% m-file must take two inputs - first, the value
of
% the independent variable; second, the value of
the
% dependent variable
% t0 initial value of the independent variable
% x0 initial value of the dependent variable(s)
% if solving a system of equations, this should
be a
% row vector containing all initial values
% tf final value of the independent variable
% N number of uniformly sized time steps to be
taken to
% advance the solution from t = t0 to t = tf
%
% output:
% wi vector / matrix containing values of the
approximate
% solution to the differential equation
% ti vector containing the values of the independent
% variable at which an approximate solution has
been
% obtained
%
neqn = length ( x0 );
ti = linspace ( t0, tf, N+1 );
wi = [ zeros( neqn, N+1 ) ];
wi(1:neqn, 1) = x0';
h = ( tf - t0 ) / N;
for i = 1:N
k1 = h * feval ( RHS, t0, x0 );
k2 = h * feval ( RHS, t0 + h/2, x0 + k1/2 );
k3 = h * feval ( RHS, t0 + h/2, x0 + k2/2 );
k4 = h * feval ( RHS, t0 + h, x0 + k3 );
x0 = x0 + ( k1 + 2*k2 + 2*k3 + k4 ) / 6;
t0 = t0 + h;
wi(1:neqn,i+1) = x0';
end;
```

Adaptive Runge-Kutta method:

```
function [wi, ti, count] = rkf45 ( RHS, t0, x0, tf, parms )
%RKF45 approximate the solution of the initial value problem
%
% x'(t) = RHS( t, x ), x(t0) = x0
%
% using Runge-Kutta-Fehlberg 4th order - 5th order method -
% this routine will work for a system of first-order
equations
% as well as for a single equation
%
% calling sequences:
% [wi, ti] = rkf45 ( RHS, t0, x0, tf, parms )
% rkf45 ( RHS, t0, x0, tf, parms )
%
% inputs:
% RHS string containing name of m-file defining the
% right-hand side of the differential equation;
the
% m-file must take two inputs - first, the value
of
% the independent variable; second, the value of
the
% dependent variable
% t0 initial value of the independent variable
% x0 initial value of the dependent variable(s)
% if solving a system of equations, this should
be a
% row vector containing all initial values
% tf final value of the independent variable
% parms three component vector of paramter values
% parm(1) minimum allowed step size
% parm(2) maximum allowed step size
% parm(3) absolute error tolerance
%
% output:
% wi vector / matrix containing values of the
approximate
% solution to the differential equation
% ti vector containing the values of the independent
% variable at which an approximate solution has
been
% obtained
% count number of function evaluations used in
advancing the
% solution from t = t0 to t = tf
%
neqn = length ( x0 );
hmin = parms(1);
hmax = parms(2);
TOL = parms(3);
ti(1) = t0;
wi(1:neqn, 1) = x0';
count = 0;
h = hmax;
i = 2;
while ( t0 < tf )
k1 = h * feval ( RHS, t0, x0 );
k2 = h * feval ( RHS, t0 + h/4, x0 + k1/4 );
k3 = h * feval ( RHS, t0 + 3*h/8, x0 + 3*k1/32 + 9*k2/32 );
k4 = h * feval ( RHS, t0 + 12*h/13, x0 + 1932*k1/2197 -
7200*k2/2197 + 7296*k3/2197 );
k5 = h * feval ( RHS, t0 + h, x0 + 439*k1/216 - 8*k2 +
3680*k3/513 - 845*k4/4104 );
k6 = h * feval ( RHS, t0 + h/2, x0 - 8*k1/27 + 2*k2 -
3544*k3/2565 + 1859*k4/4104 - 11*k5/40 );
R = max ( abs ( k1/360 - 128*k3/4275 - 2197*k4/75240 + k5/50 +
2*k6/55 ) / h );
q = 0.84 * ( TOL / R ) ^ (1/4);
count = count + 6;
if ( R < TOL )
x0 = x0 + 16*k1/135 + 6656*k3/12825 + 28561*k4/56430 - 9*k5/50
+ 2*k6/55;
% x0 = x0 + 25*k1/216 + 1408*k3/2565 + 2197*k4/4104 - k5/5;
t0 = t0 + h;
ti(i) = t0;
wi(1:neqn, i) = x0';
i = i + 1;
end;
h = min ( max ( q, 0.1 ), 4.0 ) * h;
if ( h > hmax ) h = hmax; end;
if ( t0 + h > tf )
h = tf - t0;
elseif ( h < hmin )
disp ( 'Solution requires step size smaller than minimum' );
return;
end;
end;
```

A matlab code for testing all the above mentioned methods:

```
function testRKM
x0=3;
t0=1;
tf=1.3;
h=0.05;
N=floor((tf-t0)/h);
[x4, t4]=euler('example22', t0, x0, tf, N);
xt = 6./(5-3*(0:0.01:1.25).^2);
plot((0:0.01:1.25), xt, 'k-', 'LineWidth',2);
hold on;
plot(t4,x4, 'r-x', 'LineWidth',2,...
'MarkerEdgeColor','r',...
'MarkerFaceColor','r',...
'MarkerSize',10);
hold on;
h=0.05;
N=floor((tf-t0)/h);
[x4, t4]=taylor2nd('exampletaylor22', t0, x0, tf, N);
plot(t4,x4, 'b-o', 'LineWidth',2,...
'MarkerEdgeColor','b',...
'MarkerFaceColor','b',...
'MarkerSize',10);
hold on;
h=0.05;
N=floor((tf-t0)/h);
[x4, t4]=rk4('example22', t0, x0, tf, N);
plot(t4,x4, 'g-s', 'LineWidth',2,...
'MarkerEdgeColor','g',...
'MarkerFaceColor','g',...
'MarkerSize',10);
hold on;
hmin=0.000001;
hmax = 0.5;
tol=1.0e-08;
[x45, t45, c]=rkf45('example22', t0, x0, tf, [hmin hmax tol]);
plot(t45,x45, 'm o', 'LineWidth',2,...
'MarkerEdgeColor','m',...
'MarkerFaceColor','m',...
'MarkerSize',10);
hold on;
xlim([1 1.251]);
```

Graphically results are presented in figure:

Figure: Comparison of different numerical integration methods for the time step of 0.05.

Lab 7: