Problem 1

Show[VectorPlot[{y^3-9y ,x^3-9x},{x,-5,5},{y,-5,5},VectorScale->Medium,VectorPoints->20],ListPlot[{{-0,3}},PlotStyle->Red],ListPlot[{{-3,0}},PlotStyle->Red],ListPlot[{{3,0}},PlotStyle->Red],ListPlot[{{-3,0}},PlotStyle->Red]]

Problem 2

Ω = Disk[{10,10},10];
region=RegionProduct[Ω,Line[{{0},{100}}]];
lhs=D[u[x,y,t],{t,1}]-D[u[x,y,t],{x,2}]-D[u[x,y,t],{y,2}](*-r u[x,y,t](1-u[x,y,t]/k)*);
init =u[x,y,0]==Piecewise[{{10,(x-10)^2+(y-10)^2≤ 1},{0,True}}];
Clear[result];
result=NDSolveValue[{lhs==NeumannValue[0,0≤ x≤ 20 && t>0],init},u,{x,y,t}∈ region ];
Dx[x_,y_,t_]=D[result[x,y,t],{x,1}];
Dy[x_,y_,t_]=D[result[x,y,t],{y,1}];
(Subscript[Ω, 2]) = Disk[{10,10},2];
div[x_,y_,t_]=-D[Dx[x,y,t],{x,1}]-D[Dy[x,y,t],{y,1}];
curl[x_,y_,t_]=-D[Dy[x,y,t],{x,1}]-D[Dx[x,y,t],{y,1}];
Manipulate[
Quiet[Show[RegionPlot[Ω],RegionPlot[(Subscript[Ω, 2])],
VectorPlot[{-Dx[x,y,t],-Dy[x,y,t]},{x,y} ∈Ω],PlotLabel-> {div[10,10,t],curl[10,10,t]}]]
,{t,0,100,0.01}]

Manipulate[
Quiet[Show[RegionPlot[Ω],ContourPlot[result[x,y,t],{x,y} ∈ Ω,PlotLegends->Automatic],
VectorPlot[{-Dx[x,y,t],-Dy[x,y,t]},{x,y} ∈Ω],PlotLabel-> {div[10,10,t],curl[10,10,t]}]]
,{t,0,100,0.01}]

Problem 3

e=9 10^9;
q=-1.6 10^(-19);
Q=5 10^(-19);
m=9 10^(-31);
x0=10;
y0=-10;
VectorPlot[e q Q / (Sqrt[x^2+y^2] m ){x,y},{x,-10,10} ,{y,-10,10}];
Clear[x,y];
result=NDSolve[{x''[t]==  e q Q  /(m Sqrt[x[t]^2+y[t]^2]^3) x[t],y''[t]==  e q Q/(m Sqrt[x[t]^2+y[t]^2]^3) y[t],y'[0]==0 ,x'[0]==0 ,x[0]== x0,y[0]== y0},{x,y},{t,0,5}][[1]];
x[t_]=x[t]/.result;
y[t_]=y[t]/.result;
tmax=t/.FindRoot[y[t]== 0,{t,3}];
Manipulate[Quiet[
Show[VectorPlot[e q Q / (Sqrt[x^2+y^2] m ){x,y},{x,-20,20} ,{y,-20,20}],ListPlot[{{x[t],y[t]}},PlotStyle->Red]]]
,{t,0,tmax,0.001}]

Stationary stokes flow

Ω=RegionUnion[Rectangle[{0,0},{1,1/2}],Rectangle[{1,1/10},{2,2/5}]];
plotΩ=RegionPlot[Ω,AspectRatio->Automatic]

op={(InactiveDiv[({{-μ,0},{0,-μ}}.(InactiveGrad[u[x,y], {x,y}, Inactive])), {x,y}, Inactive])+Derivative[1,0][p][x,y],(InactiveDiv[({{-μ,0},{0,-μ}}.(InactiveGrad[v[x,y], {x,y}, Inactive])), {x,y}, Inactive])+Derivative[0,1][p][x,y],Derivative[1,0][u][x,y]+Derivative[0,1][v][x,y]}/.μ->1;
pde=op=={0,0,0};
bcs={
DirichletCondition[u[x,y]==4*0.3*y*(0.5-y)/(0.41)^2,x==0.],DirichletCondition[{u[x,y]==0.,v[x,y]==0.},0<x<2],
DirichletCondition[p[x,y]==0.,x== 2]};
{xVel,yVel,pressure}=NDSolveValue[{op=={0,0,0},bcs},{u,v,p},{x,y}∈Ω,Method->{"FiniteElement","InterpolationOrder"->{u->2,v->2,p->1},"MeshOptions"->{"MaxCellMeasure"->0.0005}}];

Problem 4

plotVelocity=VectorPlot[{xVel[x,y],yVel[x,y]},{x,y}∈Ω];
Show[plotΩ,plotVelocity]

plotPressure=ContourPlot[pressure[x,y],{x,y}∈Ω];
Show[plotΩ,plotPressure]

Problem 5

DxVel[x,y]=D[xVel[x,y],x];
DyVel[x,y]=D[yVel[x,y],y];
lhs=D[u[x,y,t],t]-d D[u[x,y,t],{x,2}]-d D[u[x,y,t],{y,2}]+xVel[x,y]*D[u[x,y,t],x]+yVel[x,y]*D[u[x,y,t],y]+DxVel[x,y]*u[x,y,t] +DyVel[x,y]*u[x,y,t];
init= u[x,y,0]==Piecewise[{{100,0≤x≤0.1 && 0.3≤y≤0.5},{0,True}}];
sol=NDSolveValue[{lhs==NeumannValue[0,0≤x≤2&&0≤ y≤0.5&&t>0],init},
u,{x,y,t}∈RegionProduct[Ω,Line[{{0},{1}}]]]

Manipulate[
Quiet[DensityPlot[sol[x,y,t],{x,y}∈Ω,PlotRange->{{0,2},{0,0.5}}]],
{t,0,1,0.01}
]