stable, numerical, maclaurin series of analytic functions

[Fotios Kasolis]

If i recall correctly that was asked time ago! Given an analytic function 
f:R->R you can get an approximation of the coefficients of the maclaurin 
polynomial by a function like

function [ p ] = maclaurin (f, n)
N = min (max (2 ^ n, 1024), 1048576);
x = [ 0:N - 1 ] / N;
z = exp (2 * i * pi * x);
p = real (fft (f (z) / N));
p(find (abs (p) <= eps)) = 0;
p = p(n:-1:1);
endfunction

This will not work for singular functions. For instance, it ll return Nans if 
your function is f(x) = 1/(1-x) and garbage if there is a singularity. An 
example is

> p = maclaurin (@(x)exp (sin (x)), 10)
> x = 0:0.01:10;
> pv = polyval (p, x);
> plot (x, pv)
> semilogy (x, abs(pv-exp (sin (x))))


The algorithm is stable since FFT is stable!

Enjoy!
/Fotios