This will probably be a long answer, but stick around. It is a good question that I think many have tripped over once in their careers.
What you are seeing is not at all uncommon. In fact, it is to be expected. By the way, it is arguably not terribly appropriate to use significance tests on the coefficients, when there is no randomness in the data anyway. We use regression models in places where they are not really appropriate, and then it seems natural to apply methods based on statistical tests. But there is no randomness here.
Lets look at an example to understand what happens. I'll do it in MATLAB, using some of my own tools. Try fitting an exponential function, by sampling 101 uniform data points from 0 to 1. The nice thing about my choice of function is the series for it is known and easily written down.
X = 0:.01:1;
Y = exp(X);
There is no error at all. So the residuals in the fit will be purely lack of fit.
mdl = polyn2sympoly(polyfitn(X,Y,'constant, x'));
P = polyn2sympoly(mdl)
P =
0.874119292489012 + 1.69114271520557*x
Now the first order truncated Taylor series model would have been simply 1+x. But a first order linear regression model through those points will miss above and below the curve. In fact, it will tend to cross the curve in 2 points, and the residuals will tend to be fairly quadratic looking. Since we know the actual Taylor series, the residuals as y-yhat will be...
(1 + x + x^2/2 + x^3/6) - P
ans =
0.125880707510988 - 0.691142715205567*x + 0.5*x^2 + 0.166666666666667*x^3
So in general, the residuals will be a function that will tend to be dominated by the highest order truncated term, in this case, the quadratic term. That does not mean the constant and linear terms were estimated exactly & correctly compared to the Taylor series. If fact, they will not be. A plot of the residuals should make this clear.
Jumping ahead to a higher order model, here is a 4th order polynomial model.
mdl = polyfitn(X,Y,'constant, x,x^2,x^3,x^4');
P = polyn2sympoly(mdl)
P =
1.00004767172716 + 0.998491471851398*x + 0.510486867080189*x^2 + 0.139707694443945*x^3 + 0.0694961565648233*x^4
You can see that the lower order coefficients are now starting to look like the known values from the Taylor series, and by the way, all of those coefficients will be highly significant. I'd expect them to be so, as the Taylor series had no zero coefficients. So if I look at the ratio of the parameter coefficients to their approximate standard errors, I get unsurprisingly large numbers.
mdl.Coefficients./mdl.ParameterStd
ans =
119272.343263854 8504.47662299035 1061.5504181253 192.89918511997 193.456498197328
Again, look at a residual error plot.
It should be no surprise if this looks vaguely sinusoidal. Again, I expect the residual curve to cross the x axis in 5 points for a 4th order model. In general, n+1 points for an n'th order polynomial model. It simply MUST look vaguely sinusoidal!
Should I have done something differently? I think so. While it makes some distorted sense to fit a polynomial model to the curve, polynomials are rarely a good choice for these tasks. They tend to have numerical problems. High order models tend to go crazy if you are not careful. And low order models are simply not sufficiently accurate when you have essentially exact data.
So when is it a good idea to use a polynomial regression model? Polynomials DO have their uses. Low order polynomial models are nice if you need only encapsulate the basic behavior of your data, in a local region where the low order model is sufficiently accurate. They make the most sense when your data is so noisy that you cannot estimate more than a few terms. Too often people end up then going to a higher order model. After all, any well behaved function will have a series representation, so a polynomial model of some order SHOULD suffice! High order polynomials just have so many problems, with the fit, with accurate evaluation, with extrapolation. They even have what I call "intrapolation" problems, where we see large oscillations or ringing inside the domain of the data.
A better choice of model is often a spline of some ilk. Use an interpolating spline, or smoothing spline, or a regression spline. The exact choice should depend on your goals for this model and what you know about the function. Is there any noise at all to be smoothed out? If so, then you might need a smoothing or regression spline. Can you afford to have many knots in the spline, and can't afford any residual error at all? Then an interpolating spline is arguably right. I've had cases where even though the data was known to have noise in it, the client could not afford to modify the shape of the curve in a specific highly curved region, needing to pin down the shape of that spline carefully. In that case, an interpolating spline was indicated despite the presence of known noise.
A regression spline fit with 5 knots to the exponential yields reasonably small residual errors, however note that the residuals still have a sinusoidal shape to them. Always look at the Y axis scaling though. The limits are +/- 0.000015, so reasonably tight. The vertical bars show the equally spaced knot locations.

Again, while the residuals still show a pattern, if they had any lower order shape to them, the regression model could have removed it, thus reducing the sum of squares of the errors! And since the data had no noise in it, I cannot expect to see a noisy looking residual plot as we see in traditional regression models on noisy data.
So you saw exactly what I'd expect to see. In hindsight, it is the only way the residuals could logically look.