Subroutines in matlab

Subroutines

In this tutorial we will assume that you know how to create vectors and matrices, know how to index into them, and know about loops. For more information on those topics see one of our tutorials on vectors, matrices, vector operations, loops, plotting, or executable files.
Sometimes you want to repeat a sequence of commands, but you want to be able to do so with different vectors and matrices. One way to make this easier is through the use of subroutines. Subroutines are just like executable files, but you can pass it different vectors and matrices to use.
For example, suppose you want a subroutine to perform Gaussian elimination, and you want to be able to pass the matrix and pass the vector (This example comes from the tutorial on loops). The first line in the file has to tell matlab what variables it will pass back when and done, and what variables it needs to work with. Here we will try to find x given that Ax=b.
The routine needs the matrix A and the vector B, and it will pass back the vector x. If the name of the file is called gaussElim.m, then the first line will look like this: 
function [x] = gaussElim(A,b)
If you want to pass back more than one variable, you can include the list in the brackets with commas in between the variable names (see the second example). If you do not know how to create a file see our tutorial on executable files.

---

Here is a sample listing of the file gaussElim.m:

function [x] = gaussElim(A,b)
% File gaussElim.m
%   This subroutine will perform Gaussian elmination
%   on the matrix that you pass to it.
%   i.e., given A and b it can be used to find x,
%        Ax = b
%
%   To run this file you will need to specify several
%   things:
%   A - matrix for the left hand side.
%   b - vector for the right hand side
%
%   The routine will return the vector x.
%   ex: [x] = gaussElim(A,b)
%     this will perform Gaussian elminiation to find x.
%
%


 N = max(size(A));


% Perform Gaussian Elimination

 for j=2:N,
     for i=j:N,
        m = A(i,j-1)/A(j-1,j-1);
        A(i,:) = A(i,:) - A(j-1,:)*m;
        b(i) = b(i) - m*b(j-1);
     end
 end

% Perform back substitution

 x = zeros(N,1);
 x(N) = b(N)/A(N,N);

 for j=N-1:-1:1,
   x(j) = (b(j)-A(j,j+1:N)*x(j+1:N))/A(j,j);
 end

---

To get the vector x, you simply call the routine by name. For example, you could do the following:

>> A = [1 2 3 6; 4 3 2 3; 9 9 1 -2; 4 2 2 1]

A =

     1     2     3     6
     4     3     2     3
     9     9     1    -2
     4     2     2     1

>> b = [1 2 1 4]'

b =

     1
     2
     1
     4

>> [x] = gaussElim(A,b)


x =

    0.6809
   -0.8936
    1.8085
   -0.5532

>> 

Sometimes you want your routine to call another routine that you specify. For example, here we will demonstrate a subroutine that will approximate a D.E., y'=f(x,y), using Euler's Method. The subroutine is able to call a function, f(x,y), specified by you.
Here a subroutine is defined that will approximate a D.E. using Euler's method. If you do not know how to create a file see our tutorial on executable files.

---

Here is a sample listing of the file eulerApprox.m:

function [x,y] = eulerApprox(startx,h,endx,starty,func)
% file: eulerApprox.m
% This matlab subroutine will find the approximation to
%  a D.E. given by 
%     y' = func(x,y)
%     y(startx) = starty
%
%  To run this file you will first need to specify
%  the following:
%      startx  : the starting value for x
%      h       : the step size
%      endx    : the ending value for x
%      starty  : the initial value
%      func    : routine name to calculate the right hand 
%                side of the D.E..  This must be specified
%                as a string.
%
%   ex: [x,y] = eulerApprox(0,1,1/16,1,'f');
%       Will return the approximation of a D.E.
%       where x is from 0 to 1 in steps of 1/16.
%       The initial value is 1, and the right hand
%       side is calculated in a subroutine given by
%       f.m.
%
%  The routine will generate two vectors.  The first
%  vector is x which is the grid points starting at
%  x0=0 and have a step size h.  
%
%  The second vector is an approximation to the specified
%  D.E. 
%



x = [startx:h:endx];

y = 0*x;
y(1) = starty;

for i=2:max(size(y)),
   y(i) = y(i-1) + h*feval(func,x(i-1),y(i-1));
end

---

In this example, we will approximate the D.E. y'=1/y. To do this you will have to create a file called f.m with the following commands:

function [f] = f(x,y)
% Evaluation of right hand side of a differential
% equation.  

f = 1/y;

---

With the subroutine defined, you can call it whenever necessary. Note that when you put comments on the 2nd line, it acts as a help file. Also note that the function f.m must be specified as a string, 'f'.

>> help eulerApprox

  file: eulerApprox.m
  This matlab subroutine will find the approximation to
   a D.E. given by 
      y' = func(x,y)
      y(startx) = starty
 
   To run this file you will first need to specify
   the following:
       startx  : the starting value for x
       h       : the step size
       endx    : the ending value for x
       starty  : the initial value
       func    : routine name to calculate the right hand 
                 side of the D.E..  This must be specified
                 as a string.
 
    ex: [x,y] = eulerApprox(0,1,1/16,1,'f');
        Will return the approximation of a D.E.
        where x is from 0 to 1 in steps of 1/16.
        The initial value is 1, and the right hand
        side is calculated in a subroutine given by
        f.m.
 
   The routine will generate two vectors.  The first
   vector is x which is the grid points starting at
   x0=0 and have a step size h.  
 
   The second vector is an approximation to the specified
   D.E. 
 

>> [x,y] = eulerApprox(0,1/16,1,1,'f');
>> plot(x,y)

When the subroutine is done, it returns two vectors and stores them in x and y.