Clear[r,K,a,d,x0,y0] (* Make sure these parameters are not predefined *)

Two derivatives that together define the Lotka Volterra model:

dxt =x'[t]==r x[t](1-x[t]/K)-a x[t]y[t];
dyt =y'[t]==a x[t]y[t]-d y[t];
model={dxt,dyt,x[0]==x0, y[0]==y0}

For steady states, nullclines and vectors we need algebraic versions of these ODEs :

dx=dxt[[2]]/.{x[t]->x,y[t]->y}
dy=dyt[[2]]/.{x[t]->x,y[t]->y}

Find steady states:

steady=Solve[{dx==0,dy==0},{x,y}]

r=1;K=1;a=1;c=0.1;d=0.5; x0=1;y0=0.01;
tmax=50;

steady

vect=VectorPlot[{dx,dy},{x,0,1},{y,0,1},FrameLabel->{X,Y}]

null=ContourPlot[{dx==0,dy==0},{x,-0.01,1.05},{y,-0.01,1.05},PlotPoints->100,ContourStyle->{Red,Blue},FrameLabel->{X,Y}]

Show[null,vect,Epilog->{PointSize[0.02],Point[{x,y}/.steady]}]

nsol=NDSolve[model,{x[t],y[t]},{t,0,tmax}];
nsol[[1]]/.t->tmax

NDSolve delivers an object that can be evaluated at any t for both x and y. We pull out x[t] and y[t] :

xt[t_]=x[t]/.nsol[[1]];yt[t_]=y[t]/.nsol[[1]];

timeplt=Plot[{xt[t],yt[t]},{t,0,tmax},PlotRange->{{0,tmax},{0,1}},AxesLabel->{Time,"X,Y"},PlotLegends->{X,Y}]

traject=ParametricPlot[{xt[t],yt[t]},{t,0,tmax},AxesOrigin->{0,0},AxesLabel->{X,Y}]

Show[null,vect,traject,Epilog->{PointSize[0.02],Point[{x,y}/.steady]}]

Make a coarser vector field with equal length vectors:

vect=VectorPlot[{dx,dy},{x,0,1},{y,0,1},VectorScale->{Small,Small,None},VectorPoints->Coarse];

Show[null,vect,traject,Epilog->{PointSize[0.02],Point[{x,y}/.steady]}]

Save the numerical solution in a table called data, and provide a nice print of that table:

tpoints=Table[t,{t,0,tmax,10}];data={tpoints,xt[tpoints],yt[tpoints]};
TableForm[data,TableHeadings->{{"Time","X","Y"},None}]

Perform stability analysis  by computing the eigenvalues of the Jacobian:

jacob={{D[dx,x],D[dx,y]},{D[dy,x],D[dy,y]}}

Eigenvalues[jacob]/. steady[[1]]
Eigenvalues[jacob]/. steady[[2]]
Eigenvalues[jacob]/. steady[[3]]

Parameter continuation by plotting the non-trivial steady state a function of K:

Clear[K];steady

Plot[{steady[[1]][[1]][[2]],steady[[1]][[2]][[2]]},{K,0,2},AxesOrigin->{0,0},AxesLabel->{K,"X,Y"},PlotLabels->{X,Y}]

Fit the model to data, e.g., estimate the parameters r and a from the data made in the table printed above:

K=1;Clear[r,a];model

CostFunction[r0_?NumericQ,a0_?NumericQ]:= (    (* data[1,2,3] = t, x, y *)
Pred=NDSolve[model /.{r->r0,a-> a0},{x[t],y[t]},{t,0,Max[data[[1]]]}];
xPred=(Pred[[1]] /. t-> data[[1]])[[1]][[2]];
yPred=(Pred[[1]] /. t-> data[[1]])[[2]][[2]];
residuals =Join[data[[2]]-xPred,data[[3]]-yPred];
Total[residuals^2]
)

FindMinimum[{CostFunction[r,a],0<r<2,0<a<2},{{r,0.5},{a,0.1}}]

r=1;a=1; (* Reset free parameters *)