QuadraticTaylor[{y_, p_, t_, dt_}] := {Clip[y + p dt + 1/2 (y^2 - t) dt^2, {-4, 4}],
p + (y^2 - t) dt, t + dt, dt};

With[{Parabola = ContourPlot[x == y^2, {x, 0, 10}, {y, -5, 5},ContourStyle->Black]},
 Manipulate[T = NestList[QuadraticTaylor, {0, [Alpha], 0, 10/N}, N];
  Show[Parabola,
   ListPlot[Map[{#[[3]], #[[1]]} &, T]]], {{[Alpha], 0.92437`}, -10,
   10}, {{dt, 0.1}, 0, 1}, {{N, 10000}, 1, 1000000, 1}]]

{inf = 20, pre = 32};
pa = ParametricNDSolveValue[{y''[x] == y@x^2 - x, y[0] == 0, y'[0] == alpha},
   y, {x, 0, inf}, alpha, MaxSteps -> Infinity, WorkingPrecision -> pre];

{i = 5, init = 92437/10^5};

(alphavalue =
   NestWhile[If[pa[# + 10^(-i)][inf] < Sqrt[inf], # + 10^(-i), i++; #] &, init,
 i <= Floor[pre/2] &] // Quiet) // AbsoluteTiming
(* {36.075637, 9243754874375653/10000000000000000} *)

p1 = Plot[#, {x, 0, 20}] &@{x^(1/2), pa[alphavalue][x],
      pa[alphavalue + 10^-(pre/2)][x]} // Quiet

points = 25;
grid = Array[# &, points, {0, inf}];
difforder = 4;
(* Definition of pdetoae isn't included in this code piece,
   please find it in the link above. *)
ptoa = pdetoae[y[x], grid, difforder];
set = {y''[x] == y@x^2 - x, y[0] == 0, y[inf] == Sqrt[inf]};
ae = MapAt[#[[2 ;; -2]] &, ptoa@set, 1];

sollst =
  FindRoot[ae, {y/@ grid, Sqrt[grid]}[Transpose], WorkingPrecision -> 32][[All, -1]];
sol = ListInterpolation[sollst, grid, InterpolationOrder -> difforder];

sol'[0]
(* 1.019443775765894394411167742141 <- Not accurate at all.
   This can be improved by using a larger points and difforder,
   but not quite economic. *)
(* However, the solution i.e. value of y is accurate enough: *)
p2 = Plot[sol[x], {x, 0, inf}, PlotStyle -> {Red, Thick, Dashed}];
Show[p1, p2]

soln[y10_?NumericQ, y20_?NumericQ] :=First@NDSolve[{y1'[t]==y2[t],y2'[t]-y1[t]^2 + t == 0,
y1[0] == y10, y2[0] == y20}, {y1, y2}, {t, 0, 5}];
Plot[Evaluate[{y1[t], y2[t]} /. soln[#, #] & /@Range[-1.5, 0.4, 0.1]], {t, 0, 5},
PlotRange -> All,PlotStyle -> {Red, Green}, AxesLabel -> {"t", "y"},
PlotLegends -> {"y1", "y2"}]

ParametricPlot[Evaluate[{y1[t], y2[t]} /. soln[#, #] & /@Range[-1.5, 0.4, 0.1]], {t, 0, 3},
PlotRange -> All,PlotStyle -> {Red, Green}, AxesLabel -> {"y1", "y2"}]

restart:with(plots):
ode:=diff(y(t),t,t)=y(t)^2-t:
bcs:=y(0)=0,y(N)=sqrt(N):
N:=50:
A1:=dsolve({ode,bcs},numeric):
odeplot(A1, [[t,y(t), color=red],[t,sqrt(t), color=green,linestyle=3]],0..N,axes=boxed)

A1(0)

eq = y''[x] == y[x]^2 - x;
ibcs = {y[0] == 0, y[N1] == Sqrt[N1]};
sol = First@NDSolve[Join[{eq}, ibcs], y[x], {x, 0, N1}]
Plot[{y[x]} /. sol, {x, 0, N1}]