First we solve 10x10, 20x20, 30x30, 40x40, linear algebra problems
using the gaussian elimination program I wrote.  It's "naive" in that
it just does exactly what you would have done with no pivoting or
anything.  We also check how good the solution is by looking at A*x-b

>> n = 10; A = rand(n,n); b = rand(n,1); flops(0); x = gauss_elim(A,b); flops
ans =
        1068
>> max(abs(A*x-b))
ans =
  5.5511e-016

>> n = 20; A = rand(n,n); b = rand(n,1); flops(0); x = gauss_elim(A,b); flops
ans =
        6843

>> max(abs(A*x-b))
ans =
  2.2093e-014

>> format short e
>> n = 30; A = rand(n,n); b = rand(n,1); flops(0); x = gauss_elim(A,b); flops
ans =
       21318

>> max(abs(A*x-b))
ans =
  2.0961e-013

>> n = 40; A = rand(n,n); b = rand(n,1); flops(0); x = gauss_elim(A,b); flops
ans =
       48493

>> max(abs(A*x-b))
ans =
  1.7972e-014

Use the flop counts to find how the algorithm depends on the size of
n...
>> A = [1,10,10^2,10^3;1,20,20^2,20^3;1,30,30^2,30^3;1,40,40^2,40^3]
A =
           1          10         100        1000
           1          20         400        8000
           1          30         900       27000
           1          40        1600       64000

>> b = [1068;6843;21318;48493]
b =
        1068
        6843
       21318
       48493

Solve A*x = b:
>> A\b
ans =
 -7.0000e+000
  5.8333e+000
  3.5000e+000
  6.6667e-001

So the algorithm takes 2/3 n^3 + 7/2 n^2 + 35/6 n - 7.  (I used maple
to find that 35/6... ;-)

>> help tri_diag

  x = tri_diag(A,b) takes a tri-diagonal square matrix A and a vector b and 
  returns the solution of Ax=b.

>> clear
>> n = 10; for i=1:n, A(i,i) = rand; end, for i=1:n-1, A(i,i+1) = rand; A(i+1,i) = rand; end, b = rand(n,1);
A
>> flops(0); [x] = tri_diag(A,b); flops
ans =
   174

>> n = 20; for i=1:n, A(i,i) = rand; end, for i=1:n-1, A(i,i+1) = rand; A(i+1,i) = rand; end, b = rand(n,1);
>> flops(0); [x] = tri_diag(A,b); flops
ans =
   364

>> n = 30; for i=1:n, A(i,i) = rand; end, for i=1:n-1, A(i,i+1) = rand; A(i+1,i) = rand; end, b = rand(n,1);
>> flops(0); [x] = tri_diag(A,b); flops
ans =
   554

>> n = 40; for i=1:n, A(i,i) = rand; end, for i=1:n-1, A(i,i+1) = rand; A(i+1,i) = rand; end, b = rand(n,1);
>> flops(0); [x] = tri_diag(A,b); flops
ans =
   744

We see that the flop count is linear in n...
>> plot([10,20,30,40],[174,364,554,744],'.-')

>> A = [1,10,10^2,10^3;1,20,20^2,20^3;1,30,30^2,30^3;1,40,40^2,40^3]
>> b = [174;364;554;744];
>> A\b
ans =
   -16
    19
     0
     0

So the tridiagonal solve takes 19n - 16 operations for an nxn problem.

What is this "A\b" anyway?  It's matlab's version of doing a quick
solution to A*x=b
>> help slash

 Matrix division.
  \   Backslash or left division.
      A\B is the matrix division of A into B, which is roughly the
      same as INV(A)*B , except it is computed in a different way.
      If A is an N-by-N matrix and B is a column vector with N
      components, or a matrix with several such columns, then
      X = A\B is the solution to the equation A*X = B computed by
      Gaussian elimination. A warning message is printed if A is
      badly scaled or nearly singular.  A\EYE(SIZE(A)) produces the
      inverse of A.
      If A is an M-by-N matrix with M < or > N and B is a column
      vector with M components, or a matrix with several such columns,
      then X = A\B is the solution in the least squares sense to the
      under- or overdetermined system of equations A*X = B. The
      effective rank, K, of A is determined from the QR decomposition
      with pivoting. A solution X is computed which has at most K
      nonzero components per column. If K < N this will usually not
      be the same solution as PINV(A)*B.  A\EYE(SIZE(A)) produces a
      generalized inverse of A.
 
  /   Slash or right division.
      B/A is the matrix division of A into B, which is roughly the
      same as B*INV(A) , except it is computed in a different way.
      More precisely, B/A = (A'\B')'. See \.
 
  ./  Array right division.
      B./A denotes element-by-element division.  A and B
      must have the same dimensions unless one is a scalar.
      A scalar can be divided with anything.
 
  .\  Array left division.
      A.\B. denotes element-by-element division.  A and B
      must have the same dimensions unless one is a scalar.
      A scalar can be divided with anything.


We can compare matlab's slash to using matlab's inversion routine
to using my naive guassian elimination routine:
>> A = rand(100,100); b = rand(100,1);
>> x1 = A\b;
>> x2 = inv(A)*b;
>> x3 = gauss_elim(A,b);
We see that slash is better than inv is better than naive gaussian
elimination.  Also, it was a lot faster.
>> max(abs(A*x1-b))
ans =
  6.1617e-015
>> max(abs(A*x2-b))
ans =
  1.9817e-014
>> max(abs(A*x3-b))
ans =
  6.2550e-013

>> diary off