clear; close all
%
%  Taylor series approx for function of x
%  COEN 45 lab #5
%
%  Two cases: (1) y = sin(cx); (2) y=exp(-t/tc)*cos(omega*t)
%
n = 1;          %% order of desired series approx
%n = 2;
%n = 3;
a = 0.25;       %% expand about x=0.25
whichfun = 1;   %% sin function
%whichfun = 2;   %% decay function
tolerance = 0.01;   %% tolerance 1%
switch whichfun
    case 1;     %% use fsin funciton
        c = pi;
        fun = @fsin;
        ylab = sprintf('y=sin(%5.3fx)',c);
        tit = sprintf('%d-order Taylor series for sin',n);
        more_args = c;
        xlim = [0 1]; %% plot
        ylim = [0 1.2];   %% limit
    case 2      %% use decay function
        omega = 2*pi;
        tc = 0.2;
        fun = @decay;
        ylab = sprintf('y=exp(-t/%5.3f)*cos(%5.3f*t)',tc,omega);
        tit = sprintf('%d-order Taylor series for decay fun',n);
        more_args = [tc omega];
        xlim = [0 1];
        ylim = [-1.2 1.2];
end
%
%  Plot the function
%
xmin = xlim(1);  xmax = xlim(2);
dx = (xmax-xmin)/1000;
x = xmin:dx:xmax;
y = fun(x,0,more_args);
plot(x,y);
hold on
%
%  Find the Taylor series coeffs
%
coeff = Taylor(fun,a,n,more_args);
%
%  Set up the approx function and plot
%
yapprox = zeros(size(x));
for i=0:n
    yapprox = yapprox + coeff(i+1)*(x-a).^i;
end;
plot(x,yapprox,'r');
grid
xlabel('x');
ylabel(ylab);
legend('Exact','Approx');
plot(a,fun(a,0,more_args),'*');
title(tit);
set(gca,'XLim',xlim,'YLim',ylim)
%
%  Find error bounds and plot
%
err = abs(y-yapprox)./abs(y);
ind = find(err<tolerance);
plot(x(ind([1 1])),ylim,'k','Linewidth',2)
plot(x(ind([end end])),ylim,'k','Linewidth',2)
text(x(ind(end)),0.25,'1% error bounds')