To the first part of your question:
Just write down the definition of the variables with their conditions
t := RandomReal[5]
l := RandomReal[1]
p := RandomReal[3 t]
c := RandomReal[4 (4 + p)]
then form one list
x:={t,l,p,c}
and finally generate 10 lists
xa=Array[x&,10]
{{0.81632, 0.150935, 4.17455, 15.9769}, {4.59785, 0.0985758, 1.60776,
2.19325}, {1.96392, 0.966806, 3.78854, 26.035}, {3.16595, 0.224949, 7.58408,
3.06743}, {0.590463, 0.70165, 0.470683, 5.07829}, {2.08984, 0.402493,
0.494039, 5.56079}, {2.63064, 0.461462, 0.417904, 22.6917}, {0.0377497,
0.840235, 4.49982, 20.8807}, {4.45266, 0.48599, 3.98964, 4.57971}, {3.50745,
0.0520926, 8.12091, 0.140059}}
As to your second question:
The values of t, l, etc. are the respective rows of the transposed Matrix
txa = Transpose[xa]
{{0.81632, 4.59785, 1.96392, 3.16595, 0.590463, 2.08984, 2.63064, 0.0377497,
4.45266, 3.50745}, {0.150935, 0.0985758, 0.966806, 0.224949, 0.70165,
0.402493, 0.461462, 0.840235, 0.48599, 0.0520926}, {4.17455, 1.60776,
3.78854, 7.58408, 0.470683, 0.494039, 0.417904, 4.49982, 3.98964,
8.12091}, {15.9769, 2.19325, 26.035, 3.06743, 5.07829, 5.56079, 22.6917,
20.8807, 4.57971, 0.140059}}
For example, the first element gives the list of the values of t,
txa[[1]]
{0.81632, 4.59785, 1.96392, 3.16595, 0.590463, 2.08984, 2.63064, 0.0377497,
4.45266, 3.50745}
the second one that for l
txa[[2]]
{0.150935, 0.0985758, 0.966806, 0.224949, 0.70165, 0.402493, 0.461462, \
0.840235, 0.48599, 0.0520926}
and so on.
As a check of the formulas used in the first part you can for instance compare the theoretical value of the mean of c
Integrate[16 r + 60 r s u, {r, 0, 1}, {s, 0, 1}, {u, 0, 1}]
31/2 // N
15.5
with the experimental values (in this case 5)
Array[Mean[Array[c &, 10^3]] &, 5]
{14.7292, 15.3443, 15.7023, 16.305, 14.8933}