T=1000;
sol=NDSolve[{
x'[t] == -3(x[t]-y[t]),
y'[t] == -x[t]z[t]+26x[t]-y[t],
z'[t] == x[t]y[t]-z[t],
x[0]==4,
y[0]==4,
z[0]==20
},
{x[t], y[t],z[t]}, {t, 0,T}, MaxSteps -> (10^6)
]

ParametricPlot3D[{x[t], y[t],z[t]} /. sol, {t,0,100}, PlotPoints -> (10^4)]