I think the issue is in your NC function definition. The INPUT function call used to extract the knot values from the list is using a 3.0 format, i.e.
knot = input(scan(knots, j, ' '), 3.0);
kmax = input(scan(knots, -1, ' '), 3.0);
kmax1 = input(scan(knots, -2, ' '), 3.0);
For your GLMSELECT example where the range of the X values is larger, that format looks to work okay, but for your PHREG example where the covariates are all between 0 and 1, the 3.0 format is probably giving you knot values that are not precise enough, which throws off the evaluation of the spline basis functions, and everything breaks down from there.
To check your NC function evaluations, you can output those values and compare them to the basis function evaluations in an OUTDESIGN= data set created like in this blog post of Rick's.
... View more