%function rung_err
clear all
syms x
n = 20
f = 1/(1+x^2);
syms
df = diff(f,1)
cdf = char(df)
syms
for i = 1:n+1
    i
    df = diff(df,1)
    cdfn = char(df)
% find max (or min) of f=zero of f'
    x = fzero(cdfn,0)
% evaluate f at max
    M(i) = abs(eval(cdf))
    cdf = cdfn
end
M
z = linspace(-5,5,10000);
for n = 2:20
    h = 10/(n+1);
    x = [-5:h:5]
    c = poly(x);
    r(n+1) = max(polyval(c,z));
    r(n+1) = r(n+1)/prod([1:n+2]);
end
r
