nc=15;nr=3;cx=Table[ToExpression[StringJoin["cx",ToString[i]]],{i,1,nc}];cy=Table[ToExpression[StringJoin["cy",ToString[i]]],{i,1,nc}];rx=Table[ToExpression[StringJoin["rx",ToString[i]]],{i,1,nr}];ry=Table[ToExpression[StringJoin["ry",ToString[i]]],{i,1,nr}];coordList=Flatten[{Transpose[{cx,cy}],Transpose[{rx,ry}]}];cspeed=1;rspeed=1.1;eqns=Flatten[{Table[{{Derivative[1][cx[[i]]][t],Derivative[1][cy[[i]]][t]}==cspeed*Normalize[Sum[{rx[[j]][t]-cx[[i]][t],ry[[j]][t]-cy[[i]][t]}/((cx[[i]][t]-rx[[j]][t])^2+(cy[[i]][t]-ry[[j]][t])^2),{j,1,nr}]],cx[[i]][0]==RandomReal[{-30,30}],cy[[i]][0]==RandomReal[{-30,30}]},{i,1,nc}],Table[{{Derivative[1][rx[[i]]][t],Derivative[1][ry[[i]]][t]}==rspeed*Normalize[Sum[{rx[[i]][t]-cx[[j]][t],ry[[i]][t]-cy[[j]][t]}/((cx[[j]][t]-rx[[i]][t])^2+(cy[[j]][t]-ry[[i]][t])^2),{j,1,nc}]],rx[[i]][0]==RandomReal[{-5,5}],ry[[i]][0]==RandomReal[{-5,5}]},{i,1,nr}]}];soln=NDSolve[eqns,coordList,{t,0,200},MaxSteps->200000,PrecisionGoal->2][[1]];coordListFn[t_]:=Flatten[Table[{coordList[[i]][t],coordList[[i+1]][t]},{i,1,2*(nc+nr),2}]]cops[t_]:=Evaluate[coordListFn[t][[1;;2*nc]]/.soln[[1;;2*nc]]]robbers[t_]:=Evaluate[coordListFn[t][[2*nc+1;;2*(nc+nr)]]/.soln[[2*nc+1;;2*(nc+nr)]]]Manipulate[Show[ParametricPlot[{cops[t],robbers[t]},{t,0,tmax},PlotRange->All,PlotStyle->{Darker[Red],{Darker[Green]}}],ListPlot[{Partition[cops[0],2],Partition[robbers[0],2],Partition[cops[tmax],2],Partition[robbers[tmax],2]},PlotStyle->{Directive[PointSize[Medium],Black],Directive[PointSize[Medium],Black],Darker[Red],Darker[Green]}],Axes->False],{tmax,0.001,200,1.5}]