I used Newton interpolation for function q and Horner’s algorithm for the evaluation, here is the code which works correctly;

```
function Q(q,t)
d = q;
# create divided difference
for j = 1:length(q)-1
for i = length(q):-1:j+1
d[i] = (d[i]-d[i-1])/(η[i]-η[i-j]);
end
end
# the evaluation
p = d[end] ; dp = 0;
for i = length(q)-1:-1:1
dp = p + (t-η[i]) * dp
p = d[i] + p*(t-η[i])
end
return p, dp; #p returns to Q(t) and dp returns to Q'(t)
end
```

As an example if we take `η = [1.0, 1.3, 1.6, 1.9, 2.2]; q = [0.7651977, 0.6200860 , 0.4554022, 0.2818186, 0.1103623]`

, then `Q(q , 1.5)`

it gives me the correct result `(0.5118199942386832, -0.5578831893004117)`

.

My problem is, I would like to create two functions return to `Q(q , t)[1]`

and `Q(q , t)[2]`

like `q̃ = (q,t) -> Q(q,t)[1] ; q̃′ = (q,t) -> Q(q,t)[2]`

but the result is `(0.5118199942386832, 1.6121262983073357)`

which is wrong for `Q(q , t)[2]`

.