(* TWO-DIMENSIONS *)
Clear[F,x,y,point,xmax,xmin,ymax,ymin];

F[x_,y_]={x^2-y^2,2*x*y};  (* Enter the vector field here *)

Print["F = (",F[x,y][[1]],")i + (",F[x,y][[2]],")j"];
Print["divF (x,y)= ",D[F[x,y][[1]],x]+D[F[x,y][[2]],y]];
Print["curl F = (",D[F[x,y][[2]],x]-D[F[x,y][[1]],y],")k"];

(* THREE-DIMENSIONS *)
Clear[F,x,y,z,point,xmax,xmin,ymax,ymin];

F[x_,y_,z_]={x^2-y^2,2*x*y,z^2};  (* Enter the vector field here *)

Print["F = (",F[x,y,z][[1]],")i + (",F[x,y,z][[2]],")j + (",F[x,y,z][[3]],")j"];
Print["divF (x,y)= ",D[F[x,y,z][[1]],x]+D[F[x,y,z][[2]],y]+D[F[x,y,z][[3]],z]];
Print["curl F = (",D[F[x,y,z][[3]],y]-D[F[x,y,z][[2]],z],")i + (",D[F[x,y,z][[1]],z]-D[F[x,y,z][[3]],x],")j + (",D[F[x,y,z][[2]],x]-D[F[x,y,z][[1]],y],")k"];

(* THREE-DIMENSIONS *)
Clear[F,x,y,z,q,point,xmax,xmin,ymax,ymin];

F[x_,y_,z_]={q*x/(x^2+y^2+z^2)^(3/2),q*y/(x^2+y^2+z^2)^(3/2),q*z/(x^2+y^2+z^2)^(3/2)};  (* Enter the vector field here *)

Print["F = (",F[x,y,z][[1]],")i + (",F[x,y,z][[2]],")j + (",F[x,y,z][[3]],")j"];
Print["divF (x,y)= ",D[F[x,y,z][[1]],x]+D[F[x,y,z][[2]],y]+D[F[x,y,z][[3]],z]];
Print["curl F = (",D[F[x,y,z][[3]],y]-D[F[x,y,z][[2]],z],")i + (",D[F[x,y,z][[1]],z]-D[F[x,y,z][[3]],x],")j + (",D[F[x,y,z][[2]],x]-D[F[x,y,z][[1]],y],")k"];

Simplify[]

Needs["VectorFieldPlots`"];
Do[Clear[F,x,y,point,xmax,xmin,ymax,ymin];F[x_,y_]={x,y};point={1,-1};xmin=point[[1]]-(1/(2^n));xmax=point[[1]]+(1/(2^n));ymin=point[[2]]-(1/(2^n));ymax=point[[2]]+(1/(2^n));Clear[i];i=n;picture[i]=(Needs["VectorFieldPlots`"];VectorFieldPlots`VectorFieldPlot[F[x,y],{x,xmin,xmax},{y,ymin,ymax},Frame->True,AspectRatio->1,DisplayFunction->Identity]),{n,1,8}]
Show[picture[7],DisplayFunction->$DisplayFunction]
Print["F = ",F[x,y][[1]],"i + ",F[x,y][[2]],"j"];
Print["divF (x,y)= ",(∂_x)F[x,y][[1]]+(∂_y)F[x,y][[2]]];
Print["curl F = ",(∂_x)F[x,y][[2]]-(∂_y)F[x,y][[1]]];

Show[GraphicsGrid[{{picture[7]},{picture[8]}}]]

Needs["VectorFieldPlots`"];
Do[Clear[F,x,y,point,xmax,xmin,ymax,ymin];F[x_,y_]={y,-x};point={0,0};xmin=point[[1]]-(1/(2^n));xmax=point[[1]]+(1/(2^n));ymin=point[[2]]-(1/(2^n));ymax=point[[2]]+(1/(2^n));Clear[i];i=n;picture[i]=(Needs["VectorFieldPlots`"];VectorFieldPlots`VectorFieldPlot[F[x,y],{x,xmin,xmax},{y,ymin,ymax},Frame->True,AspectRatio->1,DisplayFunction->Identity]),{n,1,8}]
Show[picture[7],DisplayFunction->$DisplayFunction]
Print["F = ",F[x,y][[1]],"i + ",F[x,y][[2]],"j"];
Print["divF (x,y)= ",(∂_x)F[x,y][[1]]+(∂_y)F[x,y][[2]]];
Print["curl F = ",(∂_x)F[x,y][[2]]-(∂_y)F[x,y][[1]]];