An array in Mathematica is called a list.

list1 = {2, 3, 5, 7, 11}   
list2 = {a, b, c, d}   
list3 = {{a, b, c}, {d, e, f}, {g, h, i}};   
list3 // MatrixForm

Getting parts of lists

list3 
list3[[2]] 
list3[[2, 3]]

The function Table is useful for making lists

Table[Prime[n], {n, 40}]

User defined Functions. For example, potential temperature

p0 = 100000. ; κ = 287. /1004. ;

θ[p_, T_] := T ((((p/p0)))^(-κ));

θ[80000. , 
      286.69]

User defined Modules. This is analogous to a subroutine.

thing[n_]:=Module[\({list, m}, 

            list = Table[Prime[m], {m, n}]; 
Total[list]\)];

thing[10000000] // Timing

Mathematica can make contour plots.

data = 
        Table[Table[Sin[x + 3 Cos[x - y]], {x, 0, 1, 0.02}], {y, 0, 1, 
            0.02}];
ListContourPlot[data, ColorFunction -> Hue]

The function ReadList will read from a file and create a list.

data1=ReadList["/Users/ross/temp2/GATE_data", {Number, Number}]

The function ListPlot will make a plot of the data.

ListPlot[data1, Axes -> False, 
      Frame -> True, PlotJoined -> False, 
      PlotStyle -> {PointSize[0.03], Hue[0]}, AspectRatio -> 1, 
      PlotRange -> All, ImageSize -> (72) 4]

Make a new list from the old list

data2=Table[{θ[data1[[k, 1]], data1[[k, 2]]], 
          100000.  - data1[[k, 1]]}, {k, Length[data1]}]

The function Show can combine plots

Print["hello"];
\(yTicks = {{0, "1000hPa"}, {20000, "800hPa"}, {40000, "600hPa"}, {60000, "400hPa"}, {80000, "200hPa"}, {100000, "0hPa"}};\)

  
listplot1=ListPlot[data2,Axes -> False,
Frame -> True,FrameTicks -> {Automatic, yTicks, yTicks, None},FrameStyle -> {Thickness[0.002]},FrameLabel -> {"θ", ""},
PlotJoined -> False,PlotStyle -> {Hue[0, 1, 0.75], PointSize[0.03]},
PlotRange -> All,AspectRatio -> 1,ImageSize ->  (72) 4,DisplayFunction -> Identity];

  
listplot2=ListPlot[data2,Axes -> False,
Frame -> True,FrameTicks -> {Automatic, yTicks, yTicks, None},FrameStyle -> {Thickness[0.002]},FrameLabel -> {"θ", ""},
PlotJoined -> True,PlotStyle -> {Hue[0, 1, 0.25], Thickness[0.007]},
PlotRange -> All,AspectRatio -> 1,ImageSize ->  (72) 4,DisplayFunction -> Identity];

plt=Show[{listplot2, listplot1},DisplayFunction -> $DisplayFunction];

The function Export can export a graphics object

Export["/Users/ross/temp2/GATE_profile.eps", plt]

Mathematica can do symbolic operations.  
Consider (∇^2)α = β.  We can symbolically differentiate.

\(\(\(α[{x_, 
          y_}] := \((Sin[π\ \ x])\)^3\ \((\ Sin[π\ y])\)^3\ \ Cos[
          5\ Sin[x] + 7\ y]\)\(
\)
    \)\)
\(\(\(β[{x_, 
          y_}] := \((∂\_\(a, a\)α[{a, 
                  b}] + ∂\_\(b, b\)α[{a, 
                  b}])\) /. {a -> x, b -> y}\)\(
\)
    \)\)
\(Plot3D[α[{x, y}], {x, 0, 1}, {y, 0, 1}, 

      PlotLabel -> "the function α", PlotRange -> All, 
      PlotPoints -> 32]\)
\(Plot3D[β[{x, y}], {x, 0, 1}, {y, 0, 1}, 

      PlotLabel -> \*""the function β where (∇^2)α = β"", PlotRange -> All, PlotPoints -> 32]\)

Mathematica can do very nice linear algebra. 
Consider A\(TraditionalForm\`α\&~\) = \(TraditionalForm\`β\&~\) where A is a matrix.   We can write \(TraditionalForm\`α\&~\) = (A^(-1))\(TraditionalForm\`β\&~\).

n = 32; Δ = 1. /(n + 1. );

  
A = 
          Table[Table[
              If[i == j, -4, 
                If[(Abs[i - j] == 1) || (Abs[i - j] == 
                        n), 1, 0]], {i, n^2}], {j, n^2}];

A // MatrixForm;

Ainv=Inverse[N[A]];//Timing

  
γ = 
          Flatten[Table[
              Table[{x, y}, {x, Δ, 
                  1 - Δ, Δ}], {y, Δ, 1 - Δ, Δ}], 1];

βtilda = β /@ γ;

αtilda = (Δ^2) Ainv . βtilda;

  
ListPlot3D[Partition[αtilda, n], 
          PlotRange -> All]; 
  
ListPlot3D[Partition[α /@ γ - αtilda, n], 
          PlotRange -> All];

Mathematica as a programing language. 
Consider ∂(ψ/∂)t = -∂(ψ/∂)x + ω(∂^2) (ψ/∂)(x^2) Burger's equation.  Animation cell->Convert To->Quicktime

<< "/Users/ross/temp2/plotfancy.m";

im = 102;  Δx = 
        1. /(im - 2. );  Δt = 0.5 Δx;

  
ω = 0.00275;

  
x = Table[(i - 3/2) Δx, {i, im}];

  
ψ1 = 
        Table[If[(1/3) < x[[i]] < (2/3), 1, 0], {i, 
            im}]; ψ2 = ψ1;

  
Do[{
(*                           advection∂(ψ/∂)x            diffusion ω(∂^2) (ψ/∂)(x^2) *)
\(Do[ψ2[\([i]\)] = ψ1[\([i]\)] - Δt \((\(ψ1[\([i + 1]\)] - ψ1[\([i - 1]\)]\)\/\(2\ Δx\))\) + Δt\ \((ω\ \(ψ1[\([i + 1]\)] - 2\ ψ1[\([i]\)] + ψ1[\([i - 1]\)]\)\/Δx\^2)\), {i, 
                2, im - 
                  1}]; 

ψ1 = ψ2; ψ1[\([1]\)] = ψ1[\([im - 1]\)]; ψ1[\([im]\)] = ψ1[\([2]\)]; 


            If[Mod[n, 4] == 0, 
              plotfancy[Δt\ n, ψ1]]; 
            If[Total[Abs[ψ1]] > 100, Break[]];\)
},{n, 200}];

Mathematica is pretty good showing 3D polygons.

n = 12;

polys=Table[Polygon[Table[{Sin[i + j + 1], Cos[ i + j + 2], 
                    0.1 (i + j)}, {i, 3}]],{j, n}];

Show[
      Graphics3D[polys]]

∫Sin[x] \[DifferentialD]x

Det[{{a, b}, {c, d}}]

FindRoot[x == Cos[x], {x, 1}]