%plots f(x)-P_i(x) for i=0,...,N-1 and for x in [0,w]
%also plots the error bound from Taylor's theorem
%true error is in blue, Taylor's bound for error is in red

N = 4; 
w = 2;
mesh = 100;      
x = linspace(0,w,mesh);   %vector of mesh equally spaced points in [0,w]
M = zeros(N,mesh);        %Taylor polynomial values, M(i+1,j) = P_i(x(j))
M(1,:) = ones(1,mesh);    %sets P_0(x(j)) = 1 for all j

for i=2:N
    for j=1:mesh
        M(i,j) = M(i-1,j)+x(j)^(i-1)/factorial(i-1);
    end
end

E = zeros(N,mesh);        %Error values, E(i+1,j) = f(x(j))-P_i(x(j))

for i=1:N
    for j=1:mesh
        E(i,j) = f(x(j))-M(i,j);
    end
end

E_bound = zeros(N,mesh);  %Taylor's theorem upper bound on error
for i=1:N
    for j=1:mesh
        E_bound(i,j) = exp(x(j))*x(j)^i/factorial(i);
    end
end

y_max = E_bound(1,mesh);  %y-axis upper limit for figures
close all;       %closes all figure windows

for i=1:N
    subplot(1,N,i)
    plot(x,E(i,:),'b',x,E_bound(i,:),'r')  %actual vs Taylor bound error
    xlabel(['degree = ' num2str(i-1) ])
    axis([0 w 0 y_max])
    axis square
end
    



