% in the file matlec1.txt on the H drive,
% MATLAB will write everything which appears on screen
a = rand(3,3)

a =

    0.9501    0.4860    0.4565
    0.2311    0.8913    0.0185
    0.6068    0.7621    0.8214

% We can save the data in our workspace with save:
save h:\data
% Note the important distinction with diary:
% diary creates a txt file, where save creates
% a file only useful to MATLAB
clear % clear all variables in workspace
load h:\data
a

a =

    0.9501    0.4860    0.4565
    0.2311    0.8913    0.0185
    0.6068    0.7621    0.8214

b = rand(3,1) % right-hand side vector

b =

    0.4447
    0.6154
    0.7919

% a*x = b is a linear system defined by A and b
x = a\b

x =

   -0.0638
    0.6995
    0.3622

b - a*x  % residual vector to verify x is a solution

ans =

  1.0e-015 *

         0
   -0.1110
         0

% The answer is almost a zero vector, 
% up to a factor of 10^(-15)
format long e
ans

ans =

                         0
   -1.110223024625157e-016
                         0

% Internally, MATLAB works with double hardware floats,
% the default display is format short
format short e
ans

ans =

            0
 -1.1102e-016
            0

format short
ans

ans =

  1.0e-015 *

         0
   -0.1110
         0

% The second way to solve a linear system is to use
% the reduced row echelon form (rref), for this we 
% first need to stack the column b to the matrix a
ab = [a b]

ab =

    0.9501    0.4860    0.4565    0.4447
    0.2311    0.8913    0.0185    0.6154
    0.6068    0.7621    0.8214    0.7919

rref(ab)

ans =

    1.0000         0         0   -0.0638
         0    1.0000         0    0.6995
         0         0    1.0000    0.3622

% after rref on the extended matrix, we see the
% solution in the last column; let us verify by
% extracting the last column from the answer ans:
y = ans(:,4) % select all rows and the fourth column

y =

   -0.0638
    0.6995
    0.3622

b - a*y % residual

ans =

  1.0e-015 *

         0
         0
   -0.1110

% We can also take the difference with the vector x:
x - y

ans =

  1.0e-015 *

    0.0278
    0.1110
   -0.1110

% we see that y is almost equal to x
% In practice, we work with vector norms:
norm(x-y)

ans =

  1.5944e-016

% norm gives us the length of the vector
% Yet another way to solve a linear system is 
% using an LU decomposition
[l,u,p] = lu(a)

l =

    1.0000         0         0
    0.2433    1.0000         0
    0.6387    0.5843    1.0000


u =

    0.9501    0.4860    0.4565
         0    0.7731   -0.0925
         0         0    0.5839


p =

     1     0     0
     0     1     0
     0     0     1

norm(a - l*u)

ans =

  2.7756e-017

% Now we can solve a*x = b, by two backsubstitutions
z = l\b % forward substitution

z =

    0.4447
    0.5072
    0.2115

s = u\z % backward substitution

s =

   -0.0638
    0.6995
    0.3622

norm(x-s)

ans =

  1.2413e-016

diary off % closes off the file matlec1.txt