For n=4, we may take the points O(0,0), A(3,0), B(3,4), C(0,4); I don't know whether the radius (5/2 here) is minimal.
For a general n, Erdos and Anning have a construction (1945), see Erdős–Anning theorem - Wikipedia

Edit. The minimal radius is 8/sqrt(15) = 2.06...  and corresponds to a cyclic quadrilateral having the lengths of the edges: 2, 3, 2, 4.

A 3d line is determined by a point P and a direction V. 
We need a necessary and sufficient condition for the concurrency of two (non-parallel)  lines.
Let's represent P and V as lists.

restart;
conc:=(P1,V1,P2,V2) -> LinearAlgebra:-Determinant(Matrix([P1-P2,V1,V2])):
A:= [1,-2,3]:    # Your data
d1:= [1,0,-3], [2,1,-2]:
d2:= [1,0, 3], [-1,1,2]:
solve( [ conc(A, [p,q,r], d1), conc(A, [p,q,r], d2) ], [p,q,r] );
#                 [[p = 6/5*q, q = q, r = -12/5*q]]
TheLine = A, eval([p,q,r], %[]); # the wanted line
#               TheLine = [1, -2, 3], [6/5*q, q, -12/5*q]
eval(%,q=5); #Or,
#               TheLine = [1, -2, 3], [6, 5, -12]

So, you may write it as:  (x-1)/6 = (y+2)/5 = (z-3)/(-12)

- If you want just the plot:

plots:-intersectplot(f - z, g - z, x = 0 .. 1, y = 0 .. 2, z = -0.01 .. 0.01, thickness = 3)

- If you want its projection, use Ruben's approach or

plots:-implicitplot(f-g, x=0..1,y=0..2)

- For a numeric parametrization see above, or use a system of ODEs.  

Your sum can be written as

The numerator represents the finite difference of order 2*k+2 of a monic polynomial of degree 2*k+4.
So, it's a polynomial of degree 2, the coefficients depending on k. The leading coefficient is (2*k+4)!/2.
The remaining two coefficients can be computed by recurrence Δⁿ⁺¹ = Δ(Δⁿ).

You are making a mess by using two kinds of variables such as T[1] and T__1.
They are displayed almost the same, and you are usind 2D input and document mode (a good choice only for presentations).

Execute indets(equation_set);   to see the variables. Then correct your equations.

It's a bad idea to use floats in symbolic computations. The discriminant B^2-4*A*C  is very close to 0, but not 0, so the conic will not be recognized as a parabola.
Try G:=convert(G, rational)  to see this.

Hence, use rational coefficients as above; you may use evalf at the end.

Let S1 and S2 be two spheres of radii r1 and r2  and centers O1 and O2. Denote by d the distance between their centers.

If A is in the intersection of S1 and S2, the angle O1AO2 equals 

arccos((r1^2 + r2^2 - d^2)/(2*r1*r2)).      [1]

Of course this makes sense only if S1 and S2 intersects (so, abs(r1-r2) <= d <= r1 + r2).

In your case d=0.  Applying blindly [1] it results a complex number (if r1 <> r2). This is not unexpected, and I don't consider it a bug; however a warning would be nice.

Note that Maple uses some equations instead of [1].
It is strange that for r1=2, r2=3, O1=O2=[-1.-1.-1],   FindAngle gives arccos(31/12) and [1] gives arccos(13/12).
Both are complex, but are not the same. However, when O1=O2=[0,0,0] (which should not matter), FindAngle gives  arccos(13/12) too!

plot(1/x, x=0.1..0.9, tickmarks=[[seq(i/10=i/`10`,i=1..9)], default]);

DS:=proc(X::set, C::set(set)) # Dynkin system
local D:=C union {X}, S, T, ok:=false;
if not `union`(C[]) subset X then error "C must contain subsets of X" fi;
while not ok do
  ok:=true;
  for S in D do for T in D do
    if S subset T and not (T minus S) in D then D:= D union {T minus S}; ok:=false fi;
  od od
od;
D
end proc:

DS({1,2,3}, {{1},{2}});nops(%);
     {{}, {1}, {2}, {3}, {1, 2}, {1, 3}, {2, 3}, {1, 2, 3}}

                                8

DS({1,2,3,4,5}, {{1,4},{2,4}, {1}}); nops(%);
{{}, {1}, {2}, {4}, {1, 2}, {1, 4}, {2, 4}, {3, 5}, {1, 2, 4},  {1, 3, 5}, {2, 3, 5}, {3, 4, 5}, {1, 2, 3, 5}, {1, 3, 4, 5},  {2, 3, 4, 5}, {1, 2, 3, 4, 5}}

                               16

DS({1,2,3,4,5}, {{2,4,5},{1,4,5}}); nops(%);
     {{}, {1, 3}, {2, 3}, {1, 4, 5}, {2, 4, 5}, {1, 2, 3, 4, 5}}

                               6

Download geometry-parab.mw

You probably want non-collinear vectors.

restart;
P:=(p,q,r) -> [p^2-q^2-r^2 ,2*p*q, 2*p*r]:
n:=13:
pyt:=[seq](seq(seq(P(p,q,r), p=floor(sqrt(q^2+r^2))+1..n), q=1..n),r=1..n):
pyt:=select( u -> evalb(igcd(u[])=1), pyt):
for i to nops(pyt) do for j from i+1 to nops(pyt) do
 A:=pyt[i]; B:=pyt[j];
 N2:=(A[2]*B[3]-A[3]*B[2])^2+(A[1]*B[3]-A[3]*B[1])^2+(A[1]*B[2]-A[2]*B[1])^2;
 if issqr(N2) then lprint(A,B) fi
od od:

[31, 24, 12], [3, 24, 16]
[19, 48, 12], [3, 16, 24]
[19, 48, 12], [3, 96, 80]
[27, 96, 16], [13, 144, 36]
[63, 120, 20], [115, 120, 48]
[63, 120, 20], [7, 24, 60]
[107, 144, 24], [3, 16, 24]
[31, 12, 24], [3, 16, 24]
[3, 24, 16], [19, 12, 48]
[3, 24, 16], [107, 24, 144]
[131, 72, 48], [7, 24, 60]
[61, 72, 36], [83, 120, 144]
[7, 60, 24], [9, 84, 28]
[7, 60, 24], [131, 48, 72]
[7, 60, 24], [103, 96, 120]
[7, 60, 24], [63, 20, 120]
[47, 140, 40], [5, 260, 208]
[13, 144, 36], [3, 80, 96]
[11, 48, 36], [39, 64, 48]
[11, 48, 36], [119, 96, 72]
[11, 48, 36], [21, 176, 132]
[11, 48, 36], [69, 208, 156]
[39, 64, 48], [119, 96, 72]
[39, 64, 48], [21, 176, 132]
[39, 64, 48], [69, 208, 156]
[119, 96, 72], [21, 176, 132]
[119, 96, 72], [69, 208, 156]
[27, 160, 60], [3, 96, 80]
[19, 12, 48], [3, 80, 96]
[61, 36, 72], [83, 144, 120]
[11, 36, 48], [39, 48, 64]
[11, 36, 48], [119, 72, 96]
[11, 36, 48], [21, 132, 176]
[11, 36, 48], [69, 156, 208]
[39, 48, 64], [119, 72, 96]
[39, 48, 64], [21, 132, 176]
[39, 48, 64], [69, 156, 208]
[119, 72, 96], [21, 132, 176]
[119, 72, 96], [69, 156, 208]
[49, 72, 72], [83, 144, 120]
[49, 72, 72], [83, 120, 144]
[103, 120, 96], [7, 24, 60]
[103, 120, 96], [3, 80, 96]
[69, 132, 88], [33, 260, 156]
[7, 24, 60], [9, 28, 84]
[115, 48, 120], [63, 20, 120]
[103, 96, 120], [3, 96, 80]
[3, 96, 80], [13, 36, 144]
[27, 16, 96], [13, 36, 144]
[69, 88, 132], [33, 156, 260]
[3, 80, 96], [27, 60, 160]
[21, 176, 132], [69, 208, 156]
[47, 40, 140], [5, 208, 260]
[21, 132, 176], [69, 156, 208]
 

Download ALG.mw

@Zeineb Maple can help to compute it, but only if one uses the right theorem.
 

restart;
alpha := x->floor(x):
f:=x->x^2+1:
int(f(x)*diff(alpha(x),x), x=0..3);  # wrong
#                               0
f(3)*alpha(3)-f(0)*alpha(0) - int(alpha(x)*diff(f(x),x), x=0..3); # correct
#                               17

As I have mentioned, Int_[0..3] alpha d alpha   does not exisist.

If P[1], ..., P[4] are your points, consider a fifth point P[5] of unknown coordinates, say (a,b).
Then take the conic (via its quadratic equation) thru P[1], ..., P[5] and impose the condition to be a parabola.
In general there will two parabolas, but it is also possible to have 0 (!). E.g. a parabola cannot pass thru the vertices of a square.