Taylor series alternative without using derivatives The Taylor series approximates a function as a polynomial, given that you know it's value at point $a$. The coefficients of Taylor series polynomial for the function $f(x)$ are defined as:
$p[n]=\frac{f^{(n)}(a)}{n!}$
where $f^{(n)}(a)$ denotes the nth derivative of f evaluated at the point a. (The derivative of order zero of f is defined to be f itself). For example, the first 5 Taylor coefficients of $e^x$ are:
$$\langle1,\frac{1}{2},\frac{1}{6},\frac{1}{24},\frac{1}{120}\rangle$$
However, I need a power series like Taylor series that does not require computing derivatives, should work on a wide variety of functions, and must use ONLY real numbers. It can use everything else (integral, infinite sum etc.), and may not be a polynomial. I have tried Fourier series and applying polynomial regression, but it doesn't work for most functions. Laurent series is not fitting my criteria since it requires complex numbers. Are there any others?
 A: Just like the Fourier series, you can use integration to find approximate functions. The advantage here is that you can always do integration numerically so this method lends itself well to numerical methods.
Another difference from a Taylor series expansion is that it can produce the best fit curve within a range of values instead of only around one value. Of course, Taylor is guaranteed to have zero error at this value, but with the integration methods, zero error will happen at an unknown position.
Example, approximate $y(x) = \sin( \pi x)$ using a polynomial between $x=0 \ldots 1$. Lets try $$y(x) \approx c_1 x + c_3 x^3 + c_5 x^5$$

*

*First order integral $$ \int_0^1 y(x)\,{\rm d}x = \int_0^1 (c_1 x + c_3 x^3 + c_5 x^5)\,{\rm d}x $$ $$ \frac{2}{\pi} = \frac{c_1}{2} + \frac{c_3}{4} + \frac{c_5}{6} \tag{1} $$

*Second-order integral. Choose a simple function $f(x)$ that mimics $y(x)$. Use it on both sides of the integral. In our case $y(x)$ is simple enough to make $f(x)=\sin(\pi x)$ also  $$ \int_0^1 y(x) \sin(\pi x)\,{\rm d}x = \int_0^1 (c_1 x + c_3 x^3 + c_5 x^5) \sin(\pi x)\,{\rm d}x $$ $$ \frac{1}{2} = \frac{c_1}{\pi} + \left( \frac{1}{\pi} - \frac{6}{\pi^3} \right) c_3 + \left( \frac{1}{\pi} - \frac{20}{\pi^3} + \frac{120}{\pi^5} \right) c_5 \tag{2}$$

*Third-order integral. Similar to the above, but with a higher-order basis function. Three terms is all we need since we have 3 coefficients in our approximation $$ \int_0^1 y(x) \sin(2 \pi x)\,{\rm d}x = \int_0^1 (c_1 x + c_3 x^3 + c_5 x^5) \sin(2 \pi x)\,{\rm d}x $$ $$ 0 = -\frac{c_1}{2 \pi} + \left( \frac{3}{4\pi^3} - \frac{1}{2\pi} \right) c_3 - \left( \frac{1}{2\pi} - \frac{5}{2 \pi^3} + \frac{15}{4\pi^4} \right) c_5 \tag{3} $$
Take equations (1), (2), and (3) and solve for the coefficients $c_1$, $c_3$ and $c_5$.
$$ \boxed{ \sin(\pi x) \approx 3.1154891857 x - 4.8785662612 x^3 + 1.7911004686 x^5 } $$
You can plot the above in Wolfram Alpha to get

As you can see it is a pretty good approximation between 0 and 1, but also between -1 and 0 since it is an odd function. Below is a plot of the error

which is maximum near and beyond the boundary of $x=1$, and zero at four places between 0 and 1. The compromise here is that the slope at $x=0$ is different from $\pi$ so the function deviates faster near zero than $\pi x$ which is the first term of the Taylor series.
