SYMBOLIC COMPUTATION
Mathematica can easily work with algebraic formulae, for example:  (18 x + 10 x^2 + 18 y + 11 x y + y^2)^3

Expanding the expression above gives a more complicated result:

ClearAll["Global`*"]
xxyy = (((18 x+10 (x^2)+18 y+11 x y+(y^2)))^3)
Expand[xxyy]

Factoring the previous expression puts it in a much simpler form:

Factor[Expand[(((18 x+10 (x^2)+18 y+11 x y+(y^2)))^3)]]

Here is  the formula for a indefinite integral; running the cell evaluates the integral (to within a constant NOT SHOWN).
I am using the fancy dancy Palettes to make this.  You could also write:  Integrate[x^4/(x^2-1),x]--same results

∫((x^4)/((x^2)-1))\[DifferentialD]x
Integrate[x^4/(x^2-1),x]

shortcuts for:  ∫((x^4)/((x^2)-1))dx
    1   esc  intt esc
   2 ∫(Control /) d□    
   3   Use Control 6  to get (□^□) where you need it 
   4  Enter var.

∫((x^4)/((x^2) - 1))\[DifferentialD]x

∫((x^4)/((x^2) -1))\[DifferentialD]x   (* enter variables etc *)

Here is a more complicated indefinite integral:

∫(x/(c+Exp[x]))\[DifferentialD]x

Mathematica can find a closed form for some sums with 
indefinite number of members using Gosper's summation algorithm. 
For example, here the letter n is the number of members in the sum, 
and we can find:

(∑_(k=1)^n)(k^2) (1+k)
Expand[%]

We can find the explicit formula for the definite integral and then get a numerical value using N[%] 
% == Previous result

(∫_2^5)(1/((x^3)-1))\[DifferentialD]x
Simplify[%]
N[%]

We can find the explicit formula for the multiple definite integral;  Consider (∫_0^1)(∫_0^x)((x^2)+(y^2))\[DifferentialD]y\[DifferentialD]x

(∫_0^x)((x^2)+(y^2))\[DifferentialD]y
(∫_0^1)%\[DifferentialD]x
(∫_0^1)(∫_0^x)((x^2)+(y^2))\[DifferentialD]y\[DifferentialD]x

Mathematica can do the Laplace transform, for example:

LaplaceTransform[(1+2 (t^2)) Exp[-3 t],t,p]

Mathematica knows the derivatives of all standard mathematical functions:

(∂_x)ArcTan[x]
D[ArcTan[x],x]

Higher order derivatives

(∂_(x,x))ArcTan[x]
D[ArcTan[x],{x,2}]

We can do repeated derivatives (should all be the same):

(∂_({x,3}))ArcTan[x]
(∂_x)(∂_x)(∂_x)ArcTan[x]
(∂_(x,x,x))ArcTan[x]
D[ArcTan[x],{x,3}]

Mathematica uses the chain rule to simplify derivatives:

(∂_x)(2 x f[(x^2)])
(∂_x)(2 x )*f[(x^2)]+ 2 x *(∂_x)f[(x^2)]

We can find multiple derivatives:

(∂_(x1,x2))Sin[x1+2 x2]

(∂_x1)Sin[x1+2 x2]
(∂_x2)((∂_x1)Sin[x1+2 x2])

Laplace’s Equation:    (∂_(x,x))u[x,y,z] + (∂_(y,y))u[x,y,z] + (∂_(z,z))u[x,y,z] == 0

ClearAll["Global`*"]
(∂_(x,x))u[x,y,z] + (∂_(y,y))u[x,y,z] + (∂_(z,z))u[x,y,z] == 0

(∂_(x,x))u[x,y,z] + (∂_(y,y))u[x,y,z] + (∂_(z,z))u[x,y,z] == 0

One D Laplace

ClearAll["Global`*"]
DSolve[(∂_(x,x))u[x] == 0, u[x],x]
uu[x_] = u[x]/.%
uuu[x_] = %[[1]]
u[x_] = uuu[x]
u[x]

The total derivative operator is Dt[function];  if you had a function x Exp[y], the Dt of this function assumes that y (might) depend on x.
Dt[f[x]] = f’[x] Dt[x];   i.e., Dt[x] = dx.

Dt[f[x]]
Dt[x Exp[y]]

ClearAll["Global`*"]
Dt[(x^m)]
SetAttributes[m,Constant]
Dt[(x^m)]
(∂_x) (x^m)

Dt[a x^2+b,x,Constants -> {a,b}]

ClearAll["Global`*"]
SetAttributes[{a,b},Constant]
Dt[a x^2+b,x]

This constructs the sum:

(∑_(j=1)^22)((x^j)/j)

Mathematica gives the exact result for the following sum (there are no variables):

(∑_(j=1)^20)(1/(j^3))

Here is the numerical value of the previous sum value to a desired precision 
(a lot of times we only want a few decimal places).

N[%]
N[(∑_(j=1)^20)(1/(j^3)),10]

Product works just like Sum:

(∏_(j=1)^9)(x+j)
Expand[%]

This finds the limit of  (sin(x)-x)/x^3  as x approaches 0:

Limit[((Sin[x]-x)/(x^3)),x->0]

Limit[((Sin[x])/x),x->0]

Abs[-(1/6)]

You can also find approximative formulae.  This computes the power series expansion of exp(-x)*sin(2x)  
about the point x=0 up to order x^6:

ClearAll["Global`*"]
s=Series[Exp[x Sin[x]],{x,0,6}]

Convert to a normal expression (drops the Order of higher terms):

norm[x_] = Normal[s]

Plot[norm[x],{x,0,4}]

We can integrate or differentiate such power series, for example, we compute here the integral of series s:  
where by running the cell above you have defined s = Series[Exp[-x Sin[x]], {x, 0, 6}]; by just "asking for s"
you can see that it is defined; then we integrate it

s
norm
∫s\[DifferentialD]x
∫norm\[DifferentialD]x

When you manipulate a power series, it is sometimes convenient to think of the series as representing functions, 
which you can, for example, compose or invert. The variable x in cos(x) is here replaced by the power series for sin(x). 
The result is the power series for cos(sin(x)):   (we show both expansions to see if they are different:

Series[Cos[x],{x,0,6}]
Series[Cos[x],{x,0,6}]/.x->Series[Sin[x],{x,0,10}]
Series[Cos[Sin[x]],{x,0,6}]

Using %

Series[Sin[x],{x,0,10}]
Series[Cos[%],{x,0,6}]
Series[Cos[Sin[x]],{x,0,6}]

Mathematica can invert power series, using the function InverseSeries. For example, 
inverting the series for tan(x) we can obtain the series for arctan(x):

Series[Tan[x],{x,0,8}]
InverseSeries[Series[Tan[x],{x,0,8}]]  (* Series of the ArcTan  *)
Series[ArcTan[x],{x,0,8}]

We can easily work with matrix- and vector- expressions. This defines a square matrix of the 4th order:

A={{1,-10,-5,6},{-10,7,-1,8},{-5,-1,3,4},{6,8,4,-1}}

We can get this matrix in a more convenient form:

MatrixForm[A]

(({{1, -10, -5, 6}, {-10, 7, -1, 8}, {-5, -1, 3, 4}, {6, 8, 4, -1}})).(({{1, -10, -5, 6}, {-10, 7, -1, 8}, {-5, -1, 3, 4}, {6, 8, 4, -1}}))
A.A
A.A.A
MatrixPower[A,3]

Here is  (A^3 + 2 A^2 - A)    (matrix multiplication and addition)

CC =  A.A.A +2 A.A -A
MatrixForm[%]

A.A.A + 2 A.A - A

By Hand:

(({{1, -10, -5, 6}, {-10, 7, -1, 8}, {-5, -1, 3, 4}, {6, 8, 4, -1}})).(({{1, -10, -5, 6}, {-10, 7, -1, 8}, {-5, -1, 3, 4}, {6, 8, 4, -1}})).(({{1, -10, -5, 6}, {-10, 7, -1, 8}, {-5, -1, 3, 4}, {6, 8, 4, -1}})) + 2 (({{1, -10, -5, 6}, {-10, 7, -1, 8}, {-5, -1, 3, 4}, {6, 8, 4, -1}})).(({{1, -10, -5, 6}, {-10, 7, -1, 8}, {-5, -1, 3, 4}, {6, 8, 4, -1}})) - (({{1, -10, -5, 6}, {-10, 7, -1, 8}, {-5, -1, 3, 4}, {6, 8, 4, -1}}))
MatrixForm[%]

You can find the determinant of A and CC

Det[A]
Det[CC]

To inverse the matrix you can use the Inverse function:

B = Inverse[A]    (*  B = (A^(-1)) *)
MatrixForm[B=Inverse[A]]

this shows that the inverse * the orig. matrix == unit matrix:

B.A
MatrixForm[B.A]

Mathematica can compute eigenvalues and of the square matrix

N[Eigenvalues[A]]
N[Eigenvectors[A]]

ClearAll["`*"]
G = {{a,b},{c,d}};
MatrixForm[G]
eig = Eigenvalues[G]   (*  Generates Eigenvalues for G  *)
eig1 = eig[[1]]
eig2 = eig[[2]]
evg = Eigenvectors[G]  (*  Generates Eigenvectors for G  *)
evg1 = evg[[1]]
evg2 = evg[[2]]

Testing G.Eigenvector = eigenvalue * Eigenvector

rhs1 = Simplify[(G.evg1)]
lhs1 =Simplify[(eig1 *evg1)]
rhs2 = Simplify[(G.evg2)] 
lhs2 = Simplify[(eig2 *evg2)]
Expand[lhs1] == Expand[rhs1]
Expand[lhs2] == Expand[rhs2]

Solving Algebraic Equations

M can also solve algebraic equations.  Note that the equation is defined with the double equals sign ==.  This is very important.  After activating the cell, see what happens if you replace the double equals with a single equals sign.

Type in and execute:  Solve[x^3-6*x^2+11*x-6==0,x]   (*  result - an array of rules *)

sol1 = Solve[x^3-6*x^2+11*x-6==0,x]   (*  result - an array of rules *)

Need to extract each root first we remove inner set of brackets using Flatten

Enter sol2 = Flatten[sol1,1]  {NOTE:  the line above is the same as the line above the line above}

sol2 = Flatten[sol1,1]

The above is a list of rules (single set of brackets);  sol1[[1]] is the first one, etc.    So enter and evaluate: 
sol2[[1]]
sol2[[2]]
sol2[[3]]

sol2[[1]]
sol2[[2]]
sol2[[3]]

Now we extract -- we write x/.sol2[[1]]  which means to take the rule  x-> 1 and assign it’s ‘assignment’ to x (i.e., x =1)
So enter and evaluate (we want to label each solution with a different variable)
X1 = x/.sol2[[1]]
X2 = x/.sol2[[2]]
X3 = x/.sol2[[3]]

X1 = x/.sol2[[1]]
X2 = x/.sol2[[2]]
X3 = x/.sol2[[3]]
XX = {X1,X2,X3}  (* put in a list just to show you how *)

(Use Solve to apply the binomial thm) i.e., solve for the values of x such that:   a (x^2) + b x + c ==0

ClearAll["Global`*"]   (*  Cleans things up *)
sol = Solve[a (x^2) + b x + c ==0, x]

Now extract the two solutions then take their ratio and product

sol = Flatten[sol,1]   (*  Note I have repeated sol  --usually OK *)
xone = x/.sol[[1]]
xtwo = x/.sol[[2]]
ratio = xone/xtwo
product = Evaluate[xone*xtwo]

Can solve simultaneous equations:   take 1.   a x + b y == c,    2.   d x + e y == f

First - a good time to clear up our variables  Use ClearAll[“`*”]   OR ClearAll[“Global`*”]; the latter goes a little deeper.
THEN:
Enter and evaluate sol= Solve[{a x + b y == c, d x + e y == f}, {x,y}]   (*  NOTE ==   and the use of  {  } *)

ClearAll["`*"]
sol= Solve[{a x + b y == c, d x + e y == f}, {x,y}]   (*  NOTE ==   and the use of  {  } *)

extract as above (AFTER YOU FLATTEN)  NOTE  instead of % we use sol to proceed

sol = Flatten[sol,1]

x = x/.sol[[1]]
y = y/. sol[[2]]
x
y

Let’s put in some numbers:

enter:
a = 2;  b = -4; c = -10; d = -8; e = 7; f = 9;
x    (* these are the values of x and y that satisfy both of the equations: {a x + b y == c, d x + e y == f}   *)
y

a = 2;  b = -4; c = -10; d = -8; e = 7; f = 9;
x    (* these are the values of x and y that satisfy both of the equations: {a x + b y == c, d x + e y == f}   *)
y

To get decimal values for x and y use N:   N[x]   N[y]   Try it

N[x]  (*  N converts x and y into real number  *)
N[y]

Try    π and e  (to get these  hit esc  p  esc     and esc  ee  esc   OR use a pallette

N[π,100]
N[E,1000]

Test the following (equality)   DOES E = e   and does Pi  π??   Use ==

E ==E 
Pi ==π

Now try using ===

E === E
Pi === π

Test == and === with e == N[e, a big number]   do the same for ===

E == N[E,200]
E === N[E,200]

Try:  N[{x,y}];   Then let x = e, y = π

N[{x,y}]
x = E;
y = π
N[{x,y}]

Most Functions will (at least attempt to) distribute on a list,  {x,y} is a list.  Vectors are lists!

Start again      (this has changes)

a x + b y = c
d x + e y = f

form two equations for y as a function of x  -- these represent curves y[x] 

     (Solve each for y[x]  -->  y1[x] and y2[x])

ClearAll["`*"]
sol1 = Solve[a x + b y == c, y] ;  (*  NOTE ==  *)
sol2 = Solve[ d x + e y == f, y];   (*  NOTE ==  *)
sol1 = Flatten[sol1,1]
sol2 = Flatten[sol2,1]
y1[x_] = y/.sol1[[1]]
y2[x_] = y/. sol2[[1]]

sol = Solve[y1[x] == y2[x],x]
sol = Flatten[sol,1]
xintercept = x/.sol

a = 2;  b = -4; c = -10; d = -8; e = 7; f = 9;
y1[x]
y2[x]
Plot[{y1[x],y2[x]},{x,0,5},PlotLegends->"Expressions"]    (* need to sub in the numbers to plot *)

somewhat circular but we can find the intercept of these two lines:

sol3 = Solve[y1[x] == y2[x],x]
sol3[[1,1]]
Flatten[sol3]  (*  same   *)
x/.%
N[%]

Remember that intercept = ((c e-b f)/(-b d+a e))  and we still have these definitions:  a = 2;  b = -4; c = -10; d = -8; e = 7; f = 9;

xintercept
N[xintercept]

Now do the above all in one line to find x at the intercept and then find the value of y at the intercept

sol = Solve[y1[x] == y2[x],x]  (*  I've defined an intermediate "function" sol *)
interceptx = x/.sol[[1]]
intercepty1 = y1[interceptx]
intercepty2 = y2[interceptx]  (* should be the same as 1 *)
intercepty1 == intercepty2   (* redundent here but often done for logic tests *)

Go back to the beginning;   a, b, c, d, are cleared with the ClearAll

ClearAll["`*"]
sol1 = Solve[a x + b y == c, y] ;  (*  NOTE ==  *)
sol2 = Solve[ d x + e y == f, y];   (*  NOTE ==  *)
sol1 = Flatten[sol1,1]
sol2 = Flatten[sol2,1]
y1[x_] = y/.sol1[[1]]
y2[x_] = y/. sol2[[1]]

sol = Solve[y1[x] == y2[x],x];
sol = Flatten[sol,1];
xintercept = x/.sol  (*  This is the value of x for the intercept  *)

Now find the two y values at the intercept:   y1[eval. at x at intercept] and y2[eval. at x intercept]

y1[xintercept]
y2[xintercept]

Simplify[y1[xintercept]]
Simplify[y2[xintercept]]    (*  These are the values of

a = 2;  b = -4; c = -10; d = -8; e = 7; f = 9;
%%     (* This is the numerical y value at intercept  *)
xintercept  (*  This is the numerical value of x at the intercept *)

Some Graphing

Swipe Plot below  and then hit F1

Plot

To graph two or more functions on one graph, we create a list of functions.

Plot[{x Sin[2x], (x^3) - 3x + 1}, {x, -4, 4}]

Though this will be fine in some cases, options may be particularly helpful to clarify which graph is which .  To correspond to the list of functions, we create a list of styles.  Pay close attention here.

Plot[{x Sin[2x], (x^3) - 3x + 1}, {x, -4, 4}, PlotStyle ->  {Thickness[.01], Dashing[{.015}]}]

Note the curly braces.

Perhaps you prefer the more colorful look:

Plot[{x Sin[2x], (x^3) - 3x + 1}, {x, -4, 4}, PlotStyle -> {Hue[0], Hue[.6]}]

Or different looks for the two plots:

Plot[{x Sin[2x], (x^3) - 3x + 1}, {x, -4, 4}, PlotStyle -> {{Hue[.3], Thickness[.01]}, { Hue[.6], Dashing[{.02}]}}]

Executing the cell below will create all of the graphs.  Select the cell bracket which contains all of the graphics.  From the "Cell" menu, choose "Animate Selected Graphics."  Use the controls at the bottom left to change the speed and direction of the animation.  Click anywhere in the notebook to end the animation.  To add more detail, make the step size smaller, say .5.  (Edit the command so that the end will read {n,1,4,.5}.)

Do[Print[Plot[Sin[n x],{x,-2 π,2 π}]],{n,1,20,0.5`}]

Here, we'll use the command "Table" instead of "Do."  Also, to make interesting, inserted the "PlotRange" command.

Table[Plot[n Sin[x], {x, -6, 6}, PlotRange -> {-5, 5}], {n, 0, 4,0.1}]

Table[Plot[n Sin[x], {x, -6, 6}, PlotRange -> {-5, 5},PlotLabel->Style[Framed["n = "<>ToString[n]],16,Blue,Background->Lighter[Yellow]]], {n, 0, 4,0.1}]

Look at generic multiple functions plot:

Plot[{E^(.5  x),Log[x],Sin[x],ArcTan[x]},{x,0,2}]

Provide some styling to different curves:

Plot[{E^(.5  x),Log[x],Sin[x],ArcTan[x]},{x,0,2},PlotStyle->{Thickness[.01],Automatic,Red,Dashed}]

Add some more styling curve by curve

Plot[{E^(.5  x),Log[x],Sin[x],ArcTan[x]},{x,0,2},PlotStyle->{{Purple,Thick},{Red,Thin},{Green,Dashed},{Blue,DotDashed}}]

Plot[{E^(.5  x),Log[x],Sin[x],ArcTan[x]},{x,0,2},PlotStyle->{{Purple,Thick},{Red,Medium},{ColorData["HTML","DarkTurquoise"],Dashed},{Brown,DotDashed}},PlotLegends->"Expressions"]    (*  RGB[.3, .7,.2]*)

MORE GRAPHICS

 Here is a plot of sin(x^3) with x ranging from -2 to 2:

Plot[Sin[(x^3)],{x,-2,2}]

Here is more complicated plot of function Exp(-x^2) Cos(20 x):

Plot[Exp[-(x^2)] Cos[20 x],{x,-2,2}]

This plot gives the probability density for an excited state of a quantum mechanical harmonic oscillator. 
Here  HermiteH[6, x]  evaluates  the 6-th Hermite polynomial.

Plot[(((HermiteH[6,x] Exp[-((x^2)/2)]))^2),{x,-6,6}]

Mathematica chooses appropriate scales for plots, even when there are singularities:

Plot[(1/(Sin[x])),{x,0,10}]

You can give a list of functions to plot:

Plot[{BesselJ[1,x],BesselJ[2,x],BesselJ[3,x]},{x,0,15}]

You can also use Mathematica to make parametric plots.
In a parametric plot, you give both the x and y 
coordinates of each point as a function of the 
third parameter, say t. Here is the curve made 
by taking the x coordinate of each point to be sin(t)
and the y coordinate to be sin(2t)

ParametricPlot[{Sin[t],Sin[2 t]},{t,0,2 π}]

This gives a plot in polar coordinates, with  
rho(t)  as the radius function at angle θ:

rho[θ_]:=(1/2) (3 ((Cos[θ])^2)-1);
ParametricPlot[{rho[θ] Cos[θ],rho[θ] Sin[θ]},{θ,0,2 π}]

This makes a contour plot of the function  sin(x) sin(3y):

Plot3D[Sin[x] Sin[3 y],{x,-2,2},{y,-2,2}]

ContourPlot[Sin[x] Sin[3 y],{x,-2,2},{y,-2,2},Contours->20]

This gives a rectangular array of values {i,j}.     {Mod[m,n] gives the remainder on division of m by n

Mod[49,5]
Mod[{13,24,39,97,105},19]

The array is quite large, so we end the input with a semicolon to stop the  result from being printed out:

s=Table[Mod[{j,i},7],{i,20},{j,20}];

This gives a density plot of the array of values containing in s:

ListDensityPlot[s]

More Interesting:

s=Table[Exp[.1 i] Sin[j] (((i*j))^4),{i,20},{j,20}];
ListDensityPlot[s]

Go back up a couple of cells and re-evaluate ListDensityPlot [s]

Prettier:

ListDensityPlot[Table[y+Sin[x^2+y^2],{x,-4,4,0.2},{y,-4,4,0.2}],ColorFunction->Hue]

Mathematica can also make three-dimensional pictures. Here is the plot of the function sin(x y) for x and y between 0 and Pi:

Plot3D[((-1)/((x^2) + (y^2))),{x,-5,5},{y,-5,5}]

Here is 3D-picture of a cylinder:

Show[Graphics3D[Cylinder[]]]

Here is 3D-picture of a cone:

Show[Graphics3D[Cone[]]]

Try the following:  Go to: http://demonstrations.wolfram.com/     

(make a hyperlink as follows:  enter as text or code: http://demonstrations.wolfram.com/   Swipe it;

Then go to Menu Insert ; Hypelrink    (OR hit   Shift Control H),  fill in window Other notebook or URL and hit OK

Go to the url by clicking on it.

Mathematica can also plot a data list:

s=Table[{x,x Cos[x]},{x,-4 π,4 π,(π/20)}];
ListPlot[s,Joined->True]
ListPlot[s]

This makes a three-dimensional plot of an array of values:

s=Table[N[Exp[-(x^2)-(y^2)]] N[Cos[5 (Sqrt[(x^2)+(y^2)])]],{x,-2,2,0.1},{y,-2,2,0.1}];
ListPlot3D[s,PlotRange->All]

s=Exp[-(x^2)-(y^2)]Cos[5 (Sqrt[(x^2)+(y^2)])]
Plot3D[s,{x,-2,2},{y,-2,2}]    (*  put in PlotRange *)

Here is the plot of the function sin(x) sin(3y) for x and y between -2 and 2:

Plot3D[Sin[x] Sin[3 y],{x,-2,2},{y,-2,2}]

This redisplays the surface, as would it be seen from  a different view point (Play with ViewPoint  put in some - values)

Plot3D[Sin[x] Sin[3 y],{x,-2,2},{y,-2,2},ViewPoint->{1,1,1}]

Here is an example of a curve in space,  specified by a single parameter:

ParametricPlot3D[{Sin[u],Cos[u],(u/3)},{u,0,15}]

Here is a sphere, generated from a parametric representation:

ParametricPlot3D[{Sin[u] Cos[v], Sin[u] Sin[v], Cos[u]},
{u, 0, Pi}, {v, 0, 2 Pi}]

The definition and a picture of octahedron:

PolyhedronData["Octahedron"]

The definition and a picture of dodecahedron:

PolyhedronData["Dodecahedron"]

You can Stellate your favorite polyhedron    --  puts pyramids on each face.

<<PolyhedronOperations`;
Stellate[PolyhedronData["Octahedron"]]

PostFix    == //     THROW THIS SECTION AWAY!!!

Clear[x] 
(x^2 - Exp[x])//(x = 2)     (*  // says shove 2 in for x on the expression contained in (   ) and evaluate it.
                                    // is called PostFix    *)  
x   (*to show you that x now = 2; see next cell *)

Review

Clear[x]  (*  this is a precaution; otherwise, until reassigned, x stays == 2   CAUSE OF MANY ERRORS USING M  *)
x
3^100
N[%]   (*  OK, take note:   (1) N tries to give you a number for whatever you put in the (2) brackets [  ]
             in a useful format;   (3)  % stands for whatever sits in the last output line above  *)

Sqrt[2]     (*as is this holds the answer but no evaluation *)
N[Sqrt[2]]   (* this evaluates to some fixed no. of default digits *)
N[Sqrt[2],50]  (* this gives you 50 digits;    play with it in another notebook  Ctrl n gets you one; change 50 to 2000*)
√12    (* I used Esc sqrt Esc;   gives exact value   *)

Here is  the formula for a indefinite integral; running the cell evaluates the integral (to within a constant NOT SHOWN).
I am using the Palettes to make this.  You could also write:  Integrate[x^4/(x^2-1),x]--same results

∫((x^4)/((x^2)-1))\[DifferentialD]x
Integrate[x^4/(x^2-1),x]

Here is a more complicated indefinite integral:  (what's a PolyLog???)

∫(x/(c+Exp[x]))\[DifferentialD]x

Mathematica can find a closed form for some sums with 
indefinite number of members using Gosper's summation algorithm. 
For example, here the letter n is the number of members in the sum, 
and we can find:

(∑_(k=1)^n)(k^2) (1+k)

We can find the explicit formula for the definite integral and then get a numerical value using N[%] 
% == Previous result

(∫_2^5)(1/((x^3)-1))\[DifferentialD]x
N[%]

Mathematica knows the derivatives of all standard mathematical functions:

(∂_x)ArcTan[x]
D[ArcTan[x],x]

(∂_(x,x))ArcTan[x]
D[ArcTan[x],{x,2}]  (* note the 2 *)
(∂_x)((∂_x)ArcTan[x])

Define a function f.  Then create a table and display is using "TableForm" which will show the ordered pairs (x, f(x)) for 0 ≤ x ≤ 5 by steps of 0.2.

Clear["Global`*"];
f[x_] := (x^(-2));
inversesquares =Table[{x,f[x]}, {x,1,20,0.1}]
ListPlot[inversesquares,PlotRange-> {0,1}]

Use of Solve to solve three simultaneous equations for the variables x,y,z
a x + b y + c z  = d 
e x + f y +g z  = h
j x + k y + l z  = m

Solve[{a x + b y + c z == d, e x + f y +g z== h, j x + k y + l z == m}, {x, y, z}]

There are two basic methods of showing colors.  The first uses Red, Green, and Blue.  ("Log" refers to the natural logarithm.)

Plot[x Log[x+2], {x, -1.99, 3}, PlotStyle -> RGBColor[0, .4, .6], PlotRange->{-2,8}]

Try This:

Table[Graphics[{{RGBColor[.1 i,.1 j,.1 k],Disk[]},Text[Style[{.1 i,.1 j,.1 k },Black,Bold,Background->Yellow]]}],{i,0,10},{j,0,10},{k,0,10}]

A second way to show a graph in color is by using Hue[k].  One example is shown.  Try other values for k between 0 and 1, remembering your rainbow.

Plot[x Log[x+2], {x, -1.99, 3}, PlotStyle -> Hue[.75],PlotRange->{-2,8}]

Manipulate[Graphics[{Hue[n],Disk[]}],{n,0.01,0.99}]

A mathematical torus:

ParametricPlot3D[{(8+3 Cos[phi]) Cos[theta],(8+3 Cos[phi]) Sin[theta],3 Sin[phi]},{theta,0,2 π},{phi,0,2 π}]

Here is a three-dimensional plot of the Rieman zeta function zeta(z) in the complex plane. 
The height of the surface is given by  Abs(zeta(z)).  The color is related to the phase of zeta(z).

Plot3D[Abs[Zeta[x+I y]],{x,-2,6},{y,2,35},PlotPoints->30]

More 2D -- Parametrics

Primarily, one needs to enter the functions x = f(t), y = g(t) as a list of functions.  Beyond that, the coding for parametrically defined graphs should be fairly clear. Below -- inserted "AspectRatio -> Automatic" to make the x- and y-scales equal.

ParametricPlot[{Cos[t], Sin[ t]}, {t, 0, 2π},AspectRatio->Automatic]

Checking results

First we'll define a function f.

ClearAll["`*"]
f[x_]:=x^4 +12x^3 - 305x^2 - 1056x + 3220

Let's call our solution set "solns."

solns = Solve[f[x] == 0, x]

Use TableForm to display the output.  Here, this really isn't necessary because the numbers are quite nice.  Sometimes, however, the answers are much easier to read in this format.

Check these values by replacing them into the function.

solns
f[x] /. solns

x1 = x/.solns[[1]]

This shows one of the reasons that Mathematica shows the solution set in this format - to make substitution (into any equation) quite easy.  Then clear both f and solns.

Clear[f, solns]

Repeat the same steps above using "NSolve" with the function f(x) = (x^4)-6(x^3)+ 3x + 8. When you are done, clear your function and the name of your solution set.  (Alternatively, one can simply refer to the previous output rather than naming the solution set.)

Out of order?

Did you make a mistake somewhere so that "%" or even "%%"  just won't work? See if this scenario sounds familiar:  You made a spelling error in an input line with a "%".  You correct the spelling only to try to execute the corrected command on the faulty output.  When this occurs, replace the "%" with Out[ ],   putting the number of the proper output line in the square brackets.

Another method

The command "FindRoot" uses an iterative process called Newton's Method, and it requires an initial input which is fairly close to the actual root being sought.  Let's deal with the equation  x - Cos[x] == 0.

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

As another example, consider the cubic f(x) = (x^3)+5(x^2)-7. First, we'll define the function, calling it f.  (How clever!)

f[x_]:=(x^3)+5(x^2)-7
Plot[f[x],{x,-6,2}]

Plot the function.  Choose an appropriate domain so that you can see the three roots.

Since the smallest root is near x = -5, we'll use FindRoot and Newton's Method with this as an initial guess.

FindRoot[f[x] == 0, {x, -5}]

Based on your plot above, choose nice x-values as initial estimates to find the other two roots.  Note that bad choices for initial inputs (too far away, past a max or min, for example) can lead to some unexpected results.

Clear[f]

Systems

Here, we will solve a system of two equations, y = 2x - 3 and (x^2)-5x+2xy+(y^2)=6. First, we will graph the system to check for approximate solutions.  We will use a new method to do this. Before this can be done, use "Help" and load the package "ImplicitPlot."  Execute the cells to name and plot each of the graphs separately.

gr1=ImplicitPlot[(x^2)- 5x + 2x y + (y^2)== 6, {x, -4, 5}]

gr2 = Plot[2x - 3, {x, -3, 4}]

The interesting part is to combine the two graphs.  Use the command "Show."  Note that the graph shown is on the intersection of the two domains given.

Show[gr1, gr2]

To solve, we need a list of equations and a list of variables for which we wish to solve.

NSolve[{(x^2)-5x + 2x y + (y^2)==6, y == 2x -3}, {x, y}]

Plot3D[Sin[x y], {x, -3, 3}, {y, -3, 3},PlotPoints->60,ColorFunction->Hue]

Plot3D[Exp[1 - (x^2)+(y^2)], {x, -7, 7}, {y, -3, 3},PlotPoints->60,ColorFunction->Hue]

Plot3D[Cos[(x^2)+(y^2)], {x, -3, 3}, {y, -3, 3},PlotPoints->60,ColorFunction->Hue]

Plot3D[Sin[((Abs[x])/(Cos[y]))],{x,-1,2},{y,-1,1},PlotPoints->35,ColorOutput->GrayLevel]

ContourPlot[Cos[(x^2)+(y^2)], {x, -3, 3}, {y, -3, 3}, PlotPoints -> 25]

ContourPlot[Cos[(x^2)+(y^2)], {x, -3, 3}, {y, -3, 3}, PlotPoints -> 25,ContourShading->False]

DensityPlot[Cos[(x^2)+(y^2)], {x, -3, 3}, {y, -3, 3}, PlotPoints -> 50]

Plot3D[y/(x^3),  {x, -1, 1}, {y, -1, 1}]
Plot3D[y/(x^3),  {x, -1, 1}, {y, -1, 1},ViewPoint->{1.178, -1.274, 2.905}]
Plot3D[y/(x^3),  {x, -1, 1}, {y, -1, 1},ViewPoint->{2.605, -1.927, 0.975}]

On any of these graphs, it is also possible to change the point of view.You may have noticed in the options that the default "ViewPoint" is (1.3, -2.4, 2).Adding the option in the command line with alternate numbers is one way to change this, but there' s another method which is usually more fun.In any of the inputs above, place a comma inside the final closing bracket and leave the cursor there.From the "Input" menu, choose "3D ViewPoint Selector..." Put the cursor on the box, click, and drag it as desired.When you get to a position you like, release the mouse, and click on "Paste." The ViewPoint option will be entered automatically.Then execute the cell to see what you did! (Of course, if you prefer, you are welcome to use the drag bars or enter specific numbers in either spherical or cartesian form instead.)

ParametricPlot3D[{Cos[2t], Cos[t], Sin[3t]}, {t, 0 , 2π}]

If you' re confident, then the first thing to do is to load a special package.Under "Help," choose "Add-Ons." Then choose "Standard Packages," "Graphics" and "Animation." Execute the appropriate cell to load the package.After doing this, execute the following :  Can' t get this to work.

Clear["Global`*"];

SpinShow[Plot3D[Sin[x+Sin[y]], {x, -3, 3}, {y, -3, 3}], Frame -> 12]

Use Play or ListPlay to create a sound; to play--double click on the graphics

Plot[Sin[10000/t],{t,0,2}]
Play[Sin[10000/t],{t,0,2}]

Limits

Limit[((Sin[x])/x), x->0]

Limit[((2(x^3)-4x)/((x^3)+1)), x-> ∞]

Limit[1/x, x->0, Direction->1]
Limit[1/x, x->0, Direction->-1]
Plot[1/x,{x,-1,1}]
Plot[((Sin[x])/x),{x,-1,1}, PlotRange->{0,1}]

Hard way to take deriviatives and then easy way.

Limit[((Tan[x+h]-Tan[x])/h), h->0]
D[Tan[x],x]

Derivatives

Be sure to indicate the variable in each case.  (Note that this is one of the few cases where Mathematica does not require the full English word for the command.)

Compare use of D and the usual Limit definition of a derivative.  δx = Δx
after writing the code below:   esc d esc gives you δ;  esc del esc gives you ∇

ClearAll["Global`*"]
y[x_] := Exp[-x] ( (x^2) - 1);
D[y[x],x];
dy[x_] = Simplify[%];
Plot[{y[x],dy[x]},{x,-5,5}]
 Limit[((y [x+ δx]- y[x] )/(δx)),δx-> 0];
dylim[x_] = %
Simplify[dylim [x]]== dy[x]

D[x Log[x], x]
y[x_] := x Log[x]
y'[x]
D[y[x],x]

Find the next derivative and then simplify the result to see the more familiar form.

D[f[x]/g[x], x]

D[-(1/2) g (t^2)+(Subscript[v, 0]) t + (Subscript[s, 0]), t]

To find the third derivative, (((d^3)y)/(dx^3)), the coding at the end changes as follows.

D[Tan[x], {x, 3}]
z[x_] := Tan[x]
z'''[x]

Alternatively, once can define a function first and then use more standard notation.

f[x_]:=Sin[Sin[Sin[x]]]

f'[x]

f''[x]

On the palette, for the familiar derivative in single-variable calculus, use (∂_□) □.

(∂_x) ArcTan[x]

For multivariable functions, use (∂_(□,□))□.  Note that the parentheses are necessary.

(∂_(x,y))((x^3) (y^4))

f[x_]:=Sin[Sin[Sin[x]]]
(∂_x) f[x]

f[x_] = (x^(Log[x]))   (* type in x and then use ctrl 6  *)  
diffy = (∂_x)(f[x])  (* hit esc pd exc Ctrl -  NOTE no brackets  could use (  )  *)
(∂_y) f[x]

Plot[{f[x],diffy},{x,.1,2},PlotLegends->True]

Plot[{f[x],diffy},{x,.1,2},PlotRange->{{0,2},{-20,20}},PlotLegends->"Expressions"]

Implicit functions

Implicitly-defined functions require some special care.  First, we'll name our equation.

equat1 = (x^3)- 4 x y + (y^2)==12

We can't take a derivative with respect to x yet, since Mathematica would interpret y as a constant.  We must first express y as f[x]. 
NOTE--ReplaceAll:  expr /. rules applies a rule or list of rules in an attempt to transform each subpart of an expression expr; so below in equat1 it says to  replace y with f[x] (defined above) everywhere in equat1.  ANOTHER WAY to say it:  substitute in for y (everywhere) the value f[x].

equat1/.y->f[x]

Then differentiate and solve for f'[x].     NOTE %  == Out[] -- take previous OUT and put it where you have typed in %   NOTE2--Solve below did nothing

D[%,x]

Solve[%, f'[x]]

A max-min problem

Consider the function y = cos(x) over the interval (-π/2, π/2).  Find the area of the largest rectangle which may be constructed with two vertices on the x-axis and the other two vertices on the graph of the function.

First, we'll draw a picture to help understand the problem a little better, beginning with the graph of cos(x).

grcosx = Plot[Cos[x], {x, -2, 2}]

The next group of commands give something new.  We construct a typical rectangle and then show the rectangle and the graph together.

rect = {GrayLevel[.5],Rectangle[{-1.1,0},{1.1, Cos[1.1]}]};
Show[{Graphics[rect], grcosx}, Axes -> True]

We now define the area function, where the base is 2x and the height is Cos[x].

areaofrect[x_]:=2x Cos[x]

Find the derivative, and solve for the stationary point. Here, the FindRoot command will be appropriate.

FindRoot[areaofrect'[x]==0, {x, 1}]

Substitute this value of x to find the area.

areaofrect[x]/.%

rect = {GrayLevel[.5],Rectangle[{-0.860334,0},{0.860334, Cos[0.860334]}]};
Show[{Graphics[rect], grcosx}, Axes -> True]

Mathematica uses the chain rule to simplify derivatives:

(∂_x)(2 x f[(x^2)])
(∂_x)(2 x )*f[(x^2)]+ 2 x *(∂_x)f[(x^2)]

This constructs the sum:

(∑_(j=1)^22)((x^j)/j)

Mathematica gives the exact result for the following sum:

(∑_(j=1)^20)(1/(j^3))

Here is the numerical value of the previous sum value to a desired precision 
(a lot of times we only want a few decimal places).

N[(∑_(j=1)^20)(1/(j^3)),10]

Product works just like Sum:

(∏_(j=1)^9)(x+j)

This finds the limit of  (sin(x)-x)/x^3  as x approaches 0:

Limit[((Sin[x]-x)/(x^3)),x->0]

You can also find approximative formulae.  This computes the power series expansion of exp(-x)*sin(2x)  
about the point x=0 up to order x^6:

s=Series[Exp[-x Sin[x]],{x,0,6}]

Convert to a normal expression (drops the Order of higher terms):

s = Normal[s]

We can integrate or differentiate such power series, for example, we compute here the integral of series s:  
where by running the cell above you have defined s = Series[Exp[-x Sin[x]], {x, 0, 6}]; by just "asking for s"
you can see that it is defined; then we integrate it

s
∫s\[DifferentialD]x

Series[Tan[x],{x,0,5}]
InverseSeries[Series[Tan[x],{x,0,5}]]
Series[ArcTan[x],{x,0,5}]

You can use the function DSolve to find symbolic solutions of ordinary differential equations. 
Here we solve one differential equation:

DSolve[y'[x]+y[x]==1,y[x],x]

This solves two simultaneous differential equations:

DSolve[{y'[x]==-z[x],z'[x]==-y[x]},{y[x],z[x]},x]

You can add boundary conditions for differential equations 
by explicitly giving additional equations such as y[0] = 0. 
Here we solve boundary value problem  
z'[x] + y[x] = 1 - x*x,  y'[x] = z[x],  z[0] = 0,  y[1] = 0:

DSolve[{z'[x]+y[x]==1-x x,y'[x]==z[x],z[0]==0,y[1]==0},{y[x],z[x]},x]

You can evaluate the numerical value of multidimensional integrals 
(you fall back on numerical solutions when the exact solutions don't work!); consider the 2D function:
z = Sqrt[2. + x + y^2 + (x y)^3]

ClearAll["Global`*'"]
Plot3D[(Sqrt[2.+x+(y^2)+(((x y))^3)]),{x,0,1},{y,-1,1},PlotRange->{0,2.2}]
NIntegrate[(Sqrt[2.+x+(y^2)+(((x y))^3)]),{x,0,1},{y,-1,1}]

Noting that this function is pretty flat.   The integral ~ (xy Area) * (average height above x-y plane);  eyeballing the latter at ~ 1.7, and taking xy Area = 1.0 * 2.0 == 2.0, we come pretty close:

IntegralApproximate = 1.7 * 2.0

Mathematica can compute the numerical value of an infinite sum, for example:

NSum[(1/(j^3)),{j,1,∞}]

Here we find a root of a Bessel function near x=10.5:

FindRoot[BesselJ[0,x],{x,10.5}]

Below we generate a table of the numerical values  of the exponential function:

ClearAll["Global`*"]
data=Table[Exp[(x/5.)],{x,0,7}]
plotty = ListPlot[data,PlotRange->{0,5}, PlotStyle->PointSize[0.02]]

Here we find a least-squares fit to data of the form c1+c2*x+c3*x^2:

fitdata1 = Fit[data,{1,x,(x^2)},x]
fitty1 = Plot[fitdata1,{x,0,7},PlotRange->{0,5}, PlotStyle->{Red,Thickness[0.01]}]

1.09428497470253716`+0.0986337275820518044` x+0.0459481907679247836` (x^2)

Using 'Fit' we use a model function of the form:  c1 + c2*x + c3*x^3 + c4*x^5;   let's plot all of them together

fitdata2 = Fit[data,{1,x,(x^3),(x^5)},x]
fitty2 = Plot[%,{x,0,7},PlotRange->{0,5}, PlotStyle->{Purple,Thickness[0.007]}]
Show[plotty, fitty1, fitty2,Frame-> True]

FindFit is a more general fitting routine--will use another polynomial model

ClearAll["`*"]
data=Table[Exp[(x/5.)],{x,7}]
plotty = ListPlot[data,PlotRange->{0,5}, PlotStyle->PointSize[0.02]]
model = a + b x +c (x^2) + d (x^3);  (*  same model as above used in Fit *)
fit= FindFit[data,model,{a,b,c,d},x]
fitfunc = Evaluate[model/.fit]  (* extracts a,b,c,d and constructs the function that forms the fit *)
fitplot = Plot[fitfunc,{x,0,8},AxesOrigin-> {0,0}]  (* PLOTS fitfunc *)
Show[fitplot, plotty]    (* actually combines a Plot and a ListPlot!!! *)

We can fit multidimensional data, for example, 3D Consider a variable such as charge density that depends on position in 3D:   ρ[x,y,z]  (can't really plot it.  listy3D is a numerical list representing a few points (here:  10}, say, of the form{x, y, z, ρ[x,y,z]}.

We fit these values of ρ[x,y,z] with a general second order expansion in x,y,z  represented by:   

{1, x, y, z, x^2, y^2, z^2, x y, x z, y z}.

ClearAll["Global`*'"]
listy3D = {{0.1, -3.1, 2.7, 2.13}, {0.2, -3.3, 2.9, 3.57}, {0.3, -3.4, 3.8, 3.63}, {0.4, -3.5, 3.9, 4.38}, 
{0.5, -3.6, 4.7, 7.23}, {0.6, -3.7, 5.1, 5.34}, {0.7, -3.8, 8.9, 4.39}, 
{0.8, -3.9, 9.3, 8.64}, {0.9, -4.0, 5.3, 3.24}, {1.2, -4.2, 3.2, 1.27}};
listylisty = listy3D[[All, 4]]
h[x_,y_,z_] = Fit[listy3D, {1, x, y, z, x^2, y^2, z^2, x y, x z, y z}, {x, y, z}]
LL = Length[listy3D];
fitlisty={0,0,0,0,0,0,0,0,0,0};
For[i = 1, i <= LL, i++, fitlisty[[i]]= h[listy3D[[i]][[1]], listy3D[[i]][[2]],listy3D[[i]][[3]]]]
fitlisty

Now we do a 1D plot to compare the fitted values (in fitlisty) with the original data (in listy3D)

plottyfitlisty = ListPlot[fitlisty, PlotStyle->{Red,PointSize[0.02]}]
plottylisty3D = ListPlot[listylisty,PlotStyle->{Blue,PointSize[0.02]}]
Show[plottyfitlisty,plottylisty3D]

Not a bad fit.

Below gives another table of  {x, y}  pairs:

ClearAll["Global`*'"]
data=Table[{x,Exp[Sin[x]]},{x,0.,1.,0.1}]
plotty = ListPlot[data,PlotRange->{0,3}, PlotStyle->PointSize[0.02]]

Here we find a fit to this new data set generated from y = Exp[Sin[x];  model is of the form  c1+c2*sin(x)+c3*sin(2x)

model = c1+c2 Sin[x]+c3 Sin[2 x];
fit = FindFit[data,model,{c1,c2,c3},x]
fitfunc = Evaluate[model/.fit]  (* extracts a,b,c,d and constructs the function that forms the fit *)
fitplot = Plot[fitfunc,{x,0,1},AxesOrigin-> {0,0}]  (* PLOTS fitfunc *)
Show[fitplot, plotty]    (* actually combines a Plot and a ListPlot!!! *)

Below we find a numerical approximation to the minimum of the function: x*sin(x) close to  x = 2 Pi.   
The result is a list: the first entry is the minimum value of the function; the second entry
 is the point  at which the minimum occurs.

FindMinimum[x Sin[x],{x,2 π}]

You can also search for a local minimum in a function of several variables:

FindMinimum[(((1+2 x-y+3 z))^2)+(((x+y-z))^2)+(((3 x-z))^2)+(((x+y+z))^2),{x,1},{y,0},{z,-1}]

M can find with Constrained Minima; the constraints are given after the function in {   };  for example:

Minimize[1+2 x-3 y+8 z,{1<=x,x<=5,1<=y,y<=3,4<=z,z<=5,x+7 y<=10},{x,y,z}]

Break to a DE problem:

SOLVING EQUATIONS

Mathematica can find the exact solution to an arbitrary cubic equation:

Solve[(x^3)+34 x+1==0,x]

This generates the first solution to a cubic equation with symbolic parameters:

a=.;b=.;Solve[(x^3)+a (x^2)+b x+2==0,x][[1]]

sol = Solve[(x^3)+a (x^2)+b x+2==0,x]
xone = x/.sol[[1]]

If you need only the numerical solution, you can use NRoots to ask for numerical result, for example :

NRoots[3+3 x-7 (x^2)-(x^3)+2 (x^4)+3 (x^7)-3 (x^8)-(x^9)+(x^10)==0,x]

You can give Solve a list of simultaneous equations to solve.
Solve can find explicit solutions for a large class of simultaneous equations.  the a =.  line releases a 
Unset[a] does the same thing.   here is an example of Unsetting a function:

f[x_] := x^2
f[x]
f[x_] =.
f[x]

f[x_] := x^2
f[x]
ClearAll["`*"]
f[x]

Now back to Solve with a Unset (insures a is not defined)

a=.;
Solve[{(x^2)+(y^2)==1,x+y==a},{x,y}]

sol = Solve[{(x^2)+(y^2)==1,x+y==a},{x,y}]
sol1 = sol[[1]]
xx = x/.sol1
yy = y/.sol1

Here is more complicated example, when we can find only a numerical solution (Solve results in nonsense):

sol = Solve[{(x^3)+(y^3)==x y,x+y+x y==1},{x,y}]
sol1 = sol[[1]]

To find the numerical solution:

NSolve[{(x^3)+(y^3)==x y,x+y+x y==1},{x,y}]

N[Solve[{(x^3)+(y^3)==x y,x+y+x y==1},{x,y}],8]

With the help of the Reduce function you can often get full solutions,  for example:

Solve[a (x^2)+b x+c==0,x]

Reduce[a (x^2)+b x+c==0,x]

It is very important that by using Reduce, you can find out under what exact conditions a particular 
set of equations has solutions. For example:

Reduce[{a x-y+a==1,3 x+(y^2)+a==2,x-a y x+3 a==0},{x,y}]

You can use the function DSolve to find symbolic solutions 
of ordinary differential equations. 
Here we solve one differential equation:

DSolve[y'[x]+y[x]==1,y[x],x]

This solves two simultaneous differential equations:

DSolve[{y'[x]==-z[x],z'[x]==-y[x]},{y[x],z[x]},x]

You can add boundary conditions for differential equations 
by explicitly giving additional equations such as y[0] = 0. 
Here we solve boundary value problem  
z'[x] + y[x] = 1 - x*x,  y'[x] = z[x],  z[0] = 0,  y[1] = 0:

DSolve[{z'[x]+y[x]==1-x x,y'[x]==z[x],z[0]==0,y[1]==0},{y[x],z[x]},x]

You can evaluate the numerical value of multidimensional integrals 
(you fall back on numerical solutions when the exact solutions don't work!); consider the 2D function:
z = Sqrt[2. + x + y^2 + (x y)^3]

ClearAll["Global`*'"]
Plot3D[(Sqrt[2.+x+(y^2)+(((x y))^3)]),{x,0,1},{y,-1,1},PlotRange->{0,2.2}]
NIntegrate[(Sqrt[2.+x+(y^2)+(((x y))^3)]),{x,0,1},{y,-1,1}]

Noting that this function is pretty flat.   The integral ~ (xy Area) * (average height above x-y plane);  eyeballing the latter at ~ 1.7, and taking xy Area = 1.0 * 2.0 == 2.0, we come pretty close:

integralapproximate = 1.7 * 2.0

Mathematica can compute the numerical value of an infinite sum, for example:

NSum[(1/(j^3)),{j,1,∞}]

Here we find a root of a Bessel function near x=10.5:

FindRoot[BesselJ[0,x],{x,10.5}]

Below we generate a table of the numerical values  of the exponential function:

data2 = Table[{x,Exp[(x/5.)]},{x,7}]
data=Table[Exp[(x/5.)],{x,7}]
plotty = ListPlot[data,PlotRange->{0,5}, PlotStyle->PointSize[0.02]]

Fit:    Here we find a least-squares fit to data of the form c1+c2*x+c3*x^2:

fitdata1 = Fit[data,{1,x,(x^2)},x]
fitty1 = Plot[%,{x,0,7},PlotRange->{0,5}, PlotStyle->{Red,Thickness[0.01]}]

1.09428497470253716`+0.0986337275820518044` x+0.0459481907679247836` (x^2)

Using 'Fit' we find a fit of the form  c1 + c2*x + c3*x^3 + c4*x^5;   let's plot all of them together

fitdata2 = Fit[data,{1,x,(x^3),(x^5)},x]
fitty2 = Plot[%,{x,0,7},PlotRange->{0,5}, PlotStyle->{Purple,Thickness[0.01]}]
Show[plotty, fitty1, fitty2,Frame-> True]

We can fit multidimensional data, for example, 3D Consider a variable such as charge density that depends on position in 3D:   ρ[x,y,z]  (can't really plot it.  listy3D is a numerical list representing a few points (here:  10}, say, of the form{x, y, z, ρ[x,y,z]}.

We fit these values of ρ[x,y,z] with a general second order expansion in x,y,z 
 represented by:   {1, x, y, z, x^2, y^2, z^2, x y, x z, y z}.

ClearAll["Global`*'"]
listy3D = {{0.1, -3.1, 2.7, 2.13}, {0.2, -3.3, 2.9, 3.57}, {0.3, -3.4, 3.8, 3.63}, {0.4, -3.5, 3.9, 4.38}, 
{0.5, -3.6, 4.7, 7.23}, {0.6, -3.7, 5.1, 5.34}, {0.7, -3.8, 8.9, 4.39}, 
{0.8, -3.9, 9.3, 8.64}, {0.9, -4.0, 5.3, 3.24}, {1.2, -4.2, 3.2, 1.27}};
listylisty = listy3D[[All, 4]]
h[x_,y_,z_] = Fit[listy3D, {1, x, y, z, x^2, y^2, z^2, x y, x z, y z}, {x, y, z}]
LL = Length[listy3D];
fitlisty={0,0,0,0,0,0,0,0,0,0};
For[i = 1, i <= LL, i++, fitlisty[[i]]= h[listy3D[[i]][[1]], listy3D[[i]][[2]],listy3D[[i]][[3]]]]
fitlisty

Now we do a 1D plot to compare the fitted values (in fitlisty) with the original data (in listy3D)

plottyfitlisty = ListPlot[fitlisty, PlotStyle->{Red,PointSize[0.02]}]
plottylisty3D = ListPlot[listylisty,PlotStyle->{Blue,PointSize[0.02]}]
Show[plottyfitlisty,plottylisty3D]

Not a bad fit.

Below gives another table of  {x, y}  pairs:

data=Table[{x,Exp[Sin[x]]},{x,0.,1.,0.2}]

Here we find a fit to the new data set, of the form  c1+c2*sin(x)+c3*sin(2x)

Fit[data,{1,Sin[x],Sin[2 x]},x]

Below we find a numerical approximation to the minimum of the function: x*sin(x) close to  x = 2 Pi.   
The result is a list: the first entry is the minimum value of the function; the second entry
 is the point  at which the minimum occurs.

FindMinimum[x Sin[x],{x,2 π}]

You can also search for a local minimum in a function of several variables:

FindMinimum[(((1+2 x-y+3 z))^2)+(((x+y-z))^2)+(((3 x-z))^2)+(((x+y+z))^2),{x,1},{y,0},{z,-1}]

M can find with Constrained Minima; the constraints are given after the function in {   };  for example:

Minimize[1+2 x-3 y+8 z,{1<=x,x<=5,1<=y,y<=3,4<=z,z<=5,x+7 y<=10},{x,y,z}]

MORE GRAPHICS

 Here is a plot of sin(x^3) with x ranging from -2 to 2:

Plot[Sin[(x^3)],{x,-2,2}]

Here is more complicated plot of function Exp(-x^2) Cos(20 x):

Plot[Exp[-(x^2)] Cos[20 x],{x,-2,2}]

This plot gives the probability density for an excited state of a quantum mechanical harmonic oscillator. 
Here  HermiteH[6, x]  evaluates  the 6-th Hermite polynomial.

Plot[(((HermiteH[6,x] Exp[-((x^2)/2)]))^2),{x,-6,6}]

Mathematica chooses appropriate scales for plots, even when there are singularities:

Plot[(1/(Sin[x])),{x,0,10}]

You can give a list of functions to plot:

Plot[{BesselJ[1,x],BesselJ[2,x],BesselJ[3,x]},{x,0,15}]

You can also use Mathematica to make parametric plots.
In a parametric plot, you give both the x and y 
coordinates of each point as a function of the 
third parameter, say t. Here is the curve made 
by taking the x coordinate of each point to be sin(t)
and the y coordinate to be sin(2t)

Plot[{Sin[t],Sin[2t]},{t,0,2 π}]
ParametricPlot[{Sin[t],Sin[2t]},{t,0,2 π}]

Manipulate[ParametricPlot[{Sin[t],Sin[2t]},{t,0,m},PlotRange->{{-1,1},{-1,1}}],{m,.01,2 π}]
Plot[{Sin[t],Sin[2t]},{t,0,2 π}]

ParametricPlot[{{2Cos[t],2Sin[t]},{2Cos[t],Sin[t]},{Cos[t],2Sin[t]},{Cos[t],Sin[t]}},{t,0,2π}]

Manipulate[ParametricPlot[{{2Cos[t],2Sin[t]},{2Cos[t],Sin[t]},{Cos[t],2Sin[t]},{Cos[t],Sin[t]}},{t,0,k},PlotLegends->Automatic,PlotRange->{{-2,2},{-2,2}}],{k,0.01,2 π}]

This gives a plot in polar coordinates, with  
rho(t)  as the radius function at angle t:

rho[t_]:=(1/2) (3 ((Cos[t])^2)-1);
ParametricPlot[{rho[t] Cos[t],rho[t] Sin[t]},{t,0,2 π}]

This makes a contour plot of the function  sin(x) sin(3y):

ContourPlot[Sin[x] Sin[3 y],{x,-2,2},{y,-2,2},Contours->20]

Mod[n,m]     n modulo m  (remainder on division of n by m)

Mod[1,3]
Mod[2,3]
Mod[3,3]
Mod[4,3]

This gives a rectangular array of values. The array is quite 
large, so we end the input with a semicolon to stop the          
result from being printed out:

Table[Mod[j,i],{i,5},{j,4}]
s=Table[Mod[j,i],{i,20},{j,20}];

This gives a density plot of the array of 
values containing in s:

ListPlot3D[s]
ListDensityPlot[s]

Mathematica can also make three-dimensional pictures.
Here is the plot of the function sin(x y) for 
x and y between 0 and Pi:

Plot3D[Sin[x y],{x,0,π},{y,0,π}]

Here is 3D-picture of a cylinder:

Graphics3D[Cylinder[]]
Show[Graphics3D[Cylinder[]]]

Here is 3D-picture of a cone:

Show[Graphics3D[Cone[]]]

Here is 3D-picture of a Moebius strip:

ExampleData[{"Geometry3D","MoebiusStrip"}]

Here is 3D-picture of a Torus  -- M has to find it; maybe sitting in Illinois

ExampleData[{"Geometry3D","Torus"}]

Here is picture of a piece of wood  -- Again, M has to find it; maybe sitting in Illinois

ExampleData[{"Texture","Wood"}]

Mathematica can also plot a data list:

s=Table[{x,x Cos[x]},{x,-4 π,4 π,(π/20)}];
ListPlot[s,Joined->True]

This makes a three-dimensional plot of an array of values:

s=Table[N[Exp[-(x^2)-(y^2)]] N[Cos[5 (Sqrt[(x^2)+(y^2)])]],{x,-2,2,0.1},{y,-2,2,0.1}];
ListPlot3D[s,PlotRange->All]

Here is the plot of the function sin(x) sin(3y) for x and y between -2 and 2:

Plot3D[Sin[x] Sin[3 y],{x,-2,2},{y,-2,2},PlotPoints->30]

This redisplays the surface, as would it be seen from  a different view point (Play with ViewPoint)

Plot3D[Sin[x] Sin[3 y],{x,-2,2},{y,-2,2},PlotPoints->30,ViewPoint->{1.3,-2.4,1.5},AxesLabel-> {"x", "y", "z"}]

Plot3D[Sin[x] Sin[3 y],{x,-2,2},{y,-2,2},PlotPoints->30,ViewPoint->{3,-3,3},AxesLabel-> {"x", "y", "z"}]

Here is an example of a curve in space,  specified by a single parameter:

ParametricPlot3D[{Sin[u],Cos[u],(u/3)},{u,0,15}]

Here is a sphere, generated from a parametric representation:

ParametricPlot3D[{Sin[u] Cos[v], Sin[u] Sin[v], Cos[u]},
{u, 0, Pi}, {v, 0, 2 Pi}]

The definition and a picture of octahedron:

PolyhedronData["Octahedron"]

The definition and a picture of dodecahedron:

PolyhedronData["Dodecahedron"]

You can Stellate your favorite polyhedron

<<PolyhedronOperations`;
Stellate[PolyhedronData["Octahedron"]]

A mathematical torus:

ParametricPlot3D[{(8+3 Cos[phi]) Cos[theta],(8+3 Cos[phi]) Sin[theta],3 Sin[phi]},{theta,0,2 π},{phi,0,2 π}]

Here is a three-dimensional plot of the Rieman zeta function zeta(z) in the complex plane. 
The height of the surface is given by  Abs(zeta(z)).  The color is related to the phase of zeta(z).

Plot3D[Abs[Zeta[x+I y]],{x,-2,6},{y,2,35},PlotPoints->30]