xvals = -2:0.04:2;
fvals = sin(xvals) + 2*exp(-30*(xvals.^2));
fsd = std(fvals);
epssd = fsd/5;
numobs = size(xvals, 2);
eps = normrnd(0, epssd, 1, numobs);
yvals = fvals + eps;

knots1 = [-1.2, -0.4, 0.4, 1.2];
designmat1 = zeros(numobs, 4+size(knots1, 2));
designmat1(:, 1) = 1;
designmat1(:, 2) = xvals;
designmat1(:, 3) = xvals.^2;
designmat1(:, 4) = xvals.^3;
for i = 1:size(knots1, 2)
    designmat1(:, 4+i) = (subplus(xvals-knots1(i))).^3;
end

lmspline5 = regress(yvals', designmat1);
yhat5 = designmat1*lmspline5;

knots2 = [-1.8, -0.4, -0.2, 0, 0.2, 0.4, 1.8];
designmat2 = zeros(numobs, size(knots2, 2));
designmat2(:, 1:4) = designmat(:, 1:4);
for i = 1:size(knots2, 2)
    designmat2(:, 4+i) = (subplus(xvals-knots2(i))).^3;
end

lmspline7 = regress(yvals', designmat2);
yhat7 = designmat2*lmspline7;

knots3 = [-1.72, -1.48, -1.2, -0.92, -0.68, -0.4, -0.12, 0.12, 0.4, 0.68, 0.92, 1.2, 1.48, 1.72];
designmat3 = zeros(numobs, size(knots3, 2));
designmat3(:, 1:4) = designmat(:, 1:4);
for i = 1:size(knots3, 2)
    designmat3(:, 4+i) = (subplus(xvals-knots3(i))).^3;
end

lmspline15 = regress(yvals', designmat3);
yhat15 = designmat3*lmspline15;

plot(xvals, yhat7, '-k', 'LineWidth', 2)
hold on;
plot(xvals, yhat5, '-b', 'LineWidth', 2)
plot(xvals, yhat15, '-r', 'Linewidth', 2)
plot(xvals, yvals, '.k', 'MarkerSize', 16)
set(gca, 'FontSize', 22, 'XLim', [-2.1, 2.1], 'YLim', [-1.5, 2.3], 'XTick', -2:1:2, 'YTick', -1:1:2, 'TickDir', 'out', 'Box', 'off')
xlabel({'';'x'}, 'FontSize', 22)
ylabel({'y';''}, 'FontSize', 22)
set(gcf, 'Position', [100, 100, 1000, 800])