Cubic BSpline interpolation I'm trying to interpolate a curve using cubic bsplines, but I cannot get it to  exactly match the original curve. Could someone please give me some pointers. The original solution is from page #7 in here.
library(splines)
func <- function(x){return (sin(x)/(1+x*x))} # original curve
x1 <- seq(-4,4,length.out = 500) #domain is -4 to +4
f1vals <- func(x1)
k1 <- c(-7,-6,-5,-4,-3,-2,-1, 0, 1,2,3,4,5,6,7) # add 3 artificial knots at either ends
bb1 <- splineDesign(k1, x = x1, outer.ok = TRUE)

fit1 <- lm(f1vals ~ bb1)
curve(func, -4, 4)
lines(x1, fit1$fitted.values, lty = 2)

Dark line is the original curve, dotted is the interpolated. Why is it so different.

 A: It's so different because you're not using an interpolation spline, but a regression spline. You're using lm to solve the overdetermined system with 500 rows and 11 columns. If you want to follow the approach in the link you provide, you don't want to do that. You want to use an interpolation spline. This means solving a determined system, with as many equations as unknowns, and a non-zero determinant.  Here it is:
library(splines)
func <- function(x){return (sin(x) / (1 + x * x))} # original curve
# a first level of approximation
k1 <- seq(-4, 4, len = 9)
f1vals <- func(k1)
spline1 <- interpSpline(k1, f1vals)
xpred <- seq(-4, 4, len = 500) 
curve(func, -4, 4)
lines(predict(spline1, xpred), lty = 2, col = "red")
points(k1, f1vals, col = "blue") 


As you can see, the spline interpolant (red dotted line) passes through the interpolation points (blue dots). Away from the interpolation points, the agreement with the interpolated function is worse, but still pretty good. 
If we halve the spacing of the interpolation points, here's what happens:
# a second level of approximation 
k2 <- seq(-4, 4, len = 17)
f2vals <- func(k2)
spline2 <- interpSpline(k2, f2vals)
lines(predict(spline2, xpred), lty = 2, col = "green")
points(k2, f2vals, col = "purple") 


The new spline (green dotted line) is now extremely close to the original function. If you were to compute the maximum of the difference between the original function and the second spline approximation, you would see that it has became less than or equal to $\frac{1}{16}$ of the maximum difference using the first spline approximation. 
