Here is the equation for computing t-probability.
n=degrees of freedom of data

You need the gamma() function, which I gave you.

Define f(n) = 
gamma( (n+1)/2)/( sqrt(Pi*n) * gamma( n/2 ) * (1+ (x^2/n ))^( - (n+1) /2 )

Probability = area from  k to infinity = S( f(n),x , k , 40 )  
k stands in for infinity  and is adequate.
The shape is approximately the curve of error,
Some urls give a bad equation
In the 89  F5 5:t Pdf  gives the value of function
           F5 6:t Cdf  gives area between t values

These agree with the given equation above.

Total area -oo to oo is 1
========
t statistics analysis

http://helios.bto.ed.ac.uk/bto/statistics/tress4a.html 
offers an example. 
Create a folder in Var-link. In MODE, select this folder.
Select Statistics APP.  See list1   list2 ...

Fill with these values from url example.

list1  | list2 |
---------------|
500    | 470    |
460    | 460    |
510    | 460    |
470    | 480    |

F4 2:2-Var Stats

XList    list1
YList    list2
Freq     1
Category List  blank
Includ Categories {}
Enter Enter
-----

See
xbar=xmean               = 485
sigma x=sumx             = 1940
sigma (x^2)              = 942600
Sx=S Standard deviation  = 23.804761
ox=P Standard deviation  = 20.6155
n= number items          = 4

ybar                     = 467.5
sigma y                  = 1870
sigma (y^2)              = 874500
Sy                       = 9.5742
oy                       = 8.2916

=========
compute these quantities
(sigma x)^2= 1940^2 = 3763600
(sigma y)^2= 1910^2 = 3496900

(sigma x)^2/n = 3763600/4 = 940900
(sigma y)^2/n = 3648100/4 = 874225



sigma d^2x = sigma (x^2) - (sigma x)^2/n   = 942600 - 1940^2/4 = 1700
sigma d^2y = sigma (y^2) - (sigma y)^2/n   = 914300 - 1910^2/4 = 275

od^2x = sigma d^2x/(n-1) = 1700/3 = 566,.666
od^2y = sigma d^2y/(n-1) = 275/3 =  91.666
 
od2=od^2x/nx + od^2y/ny= 566.666/4 + 91.666/4= 164,582

od2 is the variance of the difference between the means.

os=sqrt(od2)=12.829

t=(xbar-ybar)/os= (485-477.5)/18.2= 1.3641

Go to F5 6:Cdf
Lower value t=1.3691   Upper value =20(any big number) Degree of freedom=

Enter Enter  Cdf  = probability difference in means  is not statistically significant.
See area =0.110744 =probability difference in means is not significant.

To compute t for P=.05  two tail Choose F5:2 Inverse  2: Inverse t
Area=0.025 =0.05/2  Deg of Freedom=6    this gives t= -2.446 for a 0.05 
two tail probability

To compute t for P=.001  Choose F5:2 Inverse  2: Inverse t
Area=.0005  Deg of Freedom=6    this gives t= -5.96 for a 0.001 two tail probability

Shade 2: Shade t   Lower Value -2.46   Upper Value=2.46    
Degree of freedom =6   Enter   See colored region = 0.95
====
Math operations have to be defined in list title region.
Assume list 2 has 4 numbers. Darken title region of list 4.
See list4={}   Type in list2^2. Enter  Now the contents are squares
of numbers found in list 2.
If you need more than 6 lists, you have to define the title as
listn before you can define a math function there, This applies to list 7 up.
If you accidently define a list title with a wrong name, delete it in var-link.

========
Program

Compute Student-t

tt()
Prgm
ClIO
{500,460.510,470}->xx    x data
{470,460,460,480}->yy    y data
dim(xx)->nx               number of x           4
dim(yy)->ny               number of y           4
sum(xx)->sxx              sum of x           1940
sum(yy)->syy              sum of x           1870
mean(xx)->mxx             mean of x           485
mean(yy)->myy             mean of y           467,5
xx*xx->xx2                xx^2              {250000
yy*yy->yy2                yy^2              {220900
sum(xx2)->sxx2            sum(xx^2)          942600
sum(yy2)->syy2            sum(yy^2)          874500
sxx^2->sxx22              sum(xx^2)^2       3763600
syy^2->syy22              sum(yy^2)^2       3496900
sxx22/nx=ssx                                 940900
syy22/ny=ssy                                 874225
sxx2-sxx^2/nx->sx2                             1700
syy2-syy^2/ny->ay2                             275      
sx2/(nx-1)->od2x                               566.666
sy2/(ny-1)->od2y                               91.666
od2x/nx+od2y/ny->od2                           164.5633           
sqrt(od2)->os                                   12,829
(mxx-myy)/os->t                                1.364096
nx+ny-2=n degree of freedom                     6
Disp "Degree f,n",n
Disp "t",t                                      1.3641            
Define h(n)=gamma((n+1)/2)/(sqrt(pi*n)*gamma(n/2))*(1-x^2/n)^-(n+1)/2)
Integral(h(n),x,t,30)->p    Probability
Disp "p",p                                    0.110744
1-p->pp
Disp "pp",pp                                   0.889256
EndPrgm