(* A Mathematica tutorial notebook created by Amol Dighe for those who are using Mathematica for the first time. Geared towards the students of Mathematical Methods course *)

2+2

(* Simple calculations *)

Simplify[ (x^3 + 5) (y^5 - Log[x])]

Expand[%]

(* Matrix algebra *)

aa = {{a,b},{c,d}}; 
MatrixForm[aa]

ff = {{f,g},{h,k}};
MatrixForm[ff]

MatrixForm[aa.ff]

MatrixForm[Inverse[aa]]

A = {{1,2,3},{4,5,6},{7,8,9}};
MatrixForm[A]

Eigenvalues[A]

Eigenvalues[A][[2]]

Eigenvectors[A]

Eigenvectors[A][[3]]

(* Derivatives, Definite / indefinite integrals *)

D[Cosh[x]^2, x]

Integrate[ Exp[-x^2], {x,0,Infinity}]

Integrate[ Exp[-x^2], {x,0,x}]

Integrate[Sin[x]^2/x, {x,0,5}]

NIntegrate[Sin[x]^2/x, {x,0,5}]

(* Plotting: 2d, 3d, contour, listplots *)

Plot[  x^4 - 10 x^2, {x,-5,5}]

Plot[ x Log[x^2], {x, -5,5}]

Plot3D[ 5 x^2 - 7x y+ 3 Exp[y], {x,0,2}, {y,0,1}]

ContourPlot[ 5 x^2 + 3 Exp[y], {x,0,2}, {y,0,1}]

(* Solving simultaneous equations *)

soln = Solve[ {x+y ==3, 5 x + 3 y^2 ==4}, {x,y}]

soln[[1]]

x/. soln[[1]]

x/. soln[[2]]

(* Solving differential equations *)

DSolve[D[D[y[x], x],x] + m omega^2 y[x] == 0, {y[x]}, {x}]

soln = DSolve[x (1-x) y''[x] + (c - (a + b + 1) x) y'[x] - a b y[x] ==0,
y[x], x]

(* Plotting solutions *)

numerical = soln/. {C[1] -> 1, C[2]-> 4, a-> 0.5, b-> 0.2, c -> 0.7}

Plot[Abs[y[x]] /.numerical, {x,0,2}]

(* Special functions *)

SphericalHarmonicY[2,1, theta, phi]

SphericalPlot3D[Abs[SphericalHarmonicY[2,2, theta, phi]], {theta,0,Pi}, {phi,0,2Pi}]

(* Hydrogen atom: radial functions *)

hydpsi[r_,n_,l_] := Exp[-  r/n]  (r/n)^l LaguerreL[n-l-1, 2l+1,r/n] ;

hydprob[r_,n_,l_] := r^2 hydpsi[r,n,l]^2;

Plot[hydpsi[r,1,0], {r,0,20}, PlotRange -> All]

Plot[hydprob[r,2,0], {r,0,20}, PlotRange  -> All]

(* Diffraction from a spherical hole *)

Plot[(BesselJ[1, x]/x)^2,{x,0,20}, PlotRange -> All]

(* Series expansions *)

Series[ Log[1+x], {x,0,2}]

Series[ Log[1+x], {x,1,2}]

(* Fitting *)

data = {{1,2},{3,E}, {5,Pi}, {6,7}};
Fit[data, {1,x, x^2},x]

f1 = a x Exp[-2 x] + b Sin[x];
fit =FindFit[data,f1, {a,b},x]

fitplot=Plot[f1 /. fit,{x,0,10}];
dataplot = ListPlot[data];
Show[{fitplot,dataplot}, PlotRange -> All]

(* Cascade radioactive decay *)

soln =DSolve[ {N1'[t] == - lambda1  * N1[t], N2'[t] == - lambda2 * N2[t] + lambda1 * N1 [t], N1[0] ==  nn, N2[0]==0}, {N1[t], N2[t]}, t]

n1[t_] := N1[t] /. soln
n2[t_] := N2[t] /. soln

nn1[t_] := n1[t] /. {nn -> 100, lambda1 ->2, lambda2 -> 1};
nn2[t_] := n2[t] /. {nn -> 100, lambda1 -> 2, lambda2 -> 1};
Plot[{nn1[t], nn2[t]}, {t, 0, 10}, PlotRange -> All]

(* Terminal velocity *)

soln = DSolve[{v'[t] == g - alpha v[t], v[0]==0}, v[t], t]

vv[t_] :=v[t] /. soln

vv[t]

Plot[{vv[t] /. {alpha->2, g -> 10}}, {t,0,10}, PlotRange -> All]

(* Predator-prey problem *)

DSolve[{r'[t] == a r[t] - alpha r[t] f[t], f'[t] == - b f[t] + beta r[t] f[t], r[0] == 200, f[0] == 10}, {r[t],f[t]}, t]

a = 1; b = 0.2; alpha = 0.02; beta = 0.001;
soln =NDSolve[{r'[t] == a r[t] - alpha r[t] f[t], f'[t] == - b f[t] + beta r[t] f[t], r[0] == 2000, f[0] == 100}, {r[t],f[t]}, {t,0,10}]

Plot[Evaluate[{r[t],f[t]} /. soln],{t,0,10}, PlotRange -> All]

(* Some QM-related examples *)

sigmax = {{0,1},{1,0}};
sigmay = {{0,-I}, {I, 0}};
sigmaz = {{1,0},{0,1}};

MatrixForm[sigmax]
MatrixForm[sigmay]
MatrixForm[sigmaz]

MatrixForm[Inverse[sigmax +  sigmay]]

ClebschGordan[{1,0},{1/2,1/2},{3/2,1/2}]

ClebschGordan[{1,0},{1/2,1/2},{1/2,1/2}]

(* Double-well 
Animation *)

LL = If[x > 0, 0, - Sqrt[2] Sin[Pi x ]];
Plot[Abs[LL]^2, {x,-1,1}]

RR =  If[x < 0, 0, Sqrt[2]Sin[Pi x ]];
Plot[Abs[RR]^2, {x,-1,1}]

S = 1/Sqrt[2] (RR + LL);
A = 1/Sqrt[2] (RR - LL);

deltaE = 1;
psi[t_] = 1/Sqrt[2] (S + Exp[- I  deltaE t] A);
Plot[Abs[psi[Pi/2]]^2, {x,-1,1}]

Animate[Plot[Abs[psi[t]]^2, {x,-1,1}], {t,0,10}]