n=2 (* tu podajemy liczbe punktów *)

x=Table["x"<>ToString[k],{k,1,n}]
(* UWAGA: w tej formie wielkości typu x1 są w postaci łańcuchów tekstowych; aby dostać wyrażenia Mathematici należy zmodyfikować kod tak:
"x"<>ToString[k]
na:
ToExpression["x"<>ToString[k]]
*)

y=Table["y"<>ToString[k],{k,1,n}]

S = Sum[1,{i,1,n}]

Sx=Sum[x[[i]],{i,1,n}]

Sy=Sum[y[[i]],{i,1,n}]

Sxx=Sum[x[[i]]^2,{i,1,n}]

Sxy=Sum[x[[i]]*y[[i]],{i,1,n}]

Syy=Sum[y[[i]]^2,{i,1,n}]

(*Alternatywny sposób, wszystkie wielkości S* otrzymujemy równocześnie *)
{
{S,Sy,Syy},
{Sx,Sxy,Sxyy},
{Sxx,Sxxy,Sxxyy}
}=
Table[

Sum[x[[i]]^p y[[i]]^q,{i,1,n}]

,{p,0,2},{q,0,2}];

Δ=S*Sxx-Sx*Sx (* Aby otrzymać symbol grecki ,,delta'' naciskamy ESC, potem duże S (Shift+S) a następnie jeszcze raz ESC
*)

a=(S*Sxy-Sx*Sy)/Δ//Simplify

b=(Sxx*Sy-Sx*Sxy)/Δ//Simplify

a t+b/.t->"x1"//Simplify (* dla dwóch punktów powinniśmy otrzymać prostą przechodzącą przez 2 te 2 punkty, co sprawdzamy tutaj *)

a t+b/.t->"x2"//Simplify