Phase Portrait for ODE with IVP

Phase portrait for any second order autonomous ODE can be found as follows.

Convert the ODE to state space. This results in 2 first order ODE's. Then call StreamPlot with these 2 equations.

Let the state variables be $x_1=x,x_2=x'(t)$, then taking derivatives w.r.t time gives $x'{_1}=x_2,x'{_2}=x''(t)=-16 x_1$. Now, using StreamPlot gives

    StreamPlot[{x2, -16 x1}, {x1, -2, 2}, {x2, -2, 2}]

To see the line that passes through the initial conditions $x_1(0)=1,x_2(0)=0.1$, add the option StreamPoints

StreamPlot[{x2, -16 x1}, {x1, -2, 2}, {x2, -5, 5}, 
 StreamPoints -> {{{{1, .1}, Red}, Automatic}}]

To verify the above is the correct phase plot, you can do

ClearAll[x, t]
ode = x''[t] + 16 x[t] == 0;
ic = {x[0] == 1, x'[0] == 1/10};
sol = x[t] /. First@(DSolve[{ode, ic}, x[t], t]);
ParametricPlot[Evaluate[{sol, D[sol, t]}], {t, 0, 3}, PlotStyle -> Red]

The advatage of phase plot, is that one does not have to solve the ODE first (so it works for nonlinear hard to solve ODE's).

All what you have to do is convert the ODE to state space and use function like StreamPlot

If you want to automate the part of converting the ODE to state space, you can also use Mathematica for that. Simply use StateSpaceModel and just read of the equations.

eq = x''[t] + 16 x[t] == 0;
ss = StateSpaceModel[{eq}, {{x[t], 0}, {x'[t], 0}}, {}, {x[t]}, t]

The above shows the A matrix in $x'=Ax$. So first row reads $x_1'(t)=x_2$ and second row reads $x'_2(t)=-16 x_1$

EquationTrekker works for me, but if you are not interested in looking at a range of solutions, it might be easier to just do it with ParametricPlot

x[t_] := -Cos[4 t]

ParametricPlot[{x[t], x'[t]} // Evaluate, {t, 0, 2 π}, 
 Axes -> False, PlotLabel -> PhaseTrajectory, Frame -> True, 
 FrameLabel -> {x[t], x'[t]}, GridLines -> Automatic]