Solve boundary value problem with NDSolve. How to print out approximations to a solution?

I solve particular boundary-value-problem for ODE with NDSolve “Shooting” method. The case is that solution is attained very slow, that seems to mean that boundary-value-problem which is supposed to be well-defined enough, in fact, is ill at some place. So i try figure out. First step is to see concrete values of a produced approximations to a solution. What language constructs should i use for that?

Simple example. Suppose we consider particle motion in vertical plane. We throw our particle from initial point with coordinates {x0,y0} := {0,0} and initial trajectory angle 45 degrees. And try to achieve point with coordinates {x1,y1} := {1,0} by varying particle initial velocity. We don’t know two things here: initial velocity and a duration of a motion. Here is how that toy problem can be presented and solved in mathematica:

gravity = 10; bvpsol = NDSolve[     {      {       (* ODE system (5th order) *)       {x''[u] / c[u]^2 == 0,        y''[u] / c[u]^2 == -gravity,        c'[u] == 0},       (* boundary conditions (5 items) *)       {x[0] == y[0] == 0,        x[1] == 1,        y[1] == 0,        x'[0] == y'[0]}       }      }, {x[u], y[u], c[u]}, u,      Method -> {"Shooting",        "StartingInitialConditions" -> {x[0] == y[0] == 0,          x'[0] == y'[0] == 1, c[0] == 1}}] // Flatten;  {dxdu, dydu} = D[{x[u], y[u]} /. bvpsol, u]; {vx0, vy0} = ({dxdu, dydu} / c[u] /. bvpsol) /. {u -> 0}; duration = c[u] /. bvpsol /. {u -> RandomReal[]};  ivpsol = NDSolve[{     (* ODE system *)     {x''[t] == 0, y''[t] == -gravity},     (* initial values *)     {x[0] == y[0] == 0, x'[0] == vx0, y'[0] == vy0}     }, {x[t], y[t]}, {t, 0, duration}];  ParametricPlot[{x[t], y[t]} /. ivpsol, {t, 0, duration},   GridLines -> Automatic, AspectRatio -> 1/2] 

comprehensive parabolic trajectory

Question: Now what options or language construct should i use to see approximations which are produced NDSolve while solving boundary-value-problem?