function amath3132hw4

  clear all; close all;

  fprintf('Ans. to Q(1)\n')
  disp('hit ENTER to proceed')
  pause
  
  f1=inline('x+3*cos(x)-exp(x)');
  df1=inline('1-3*sin(x)-exp(x)');
  [n,x_bisec,fx]=bisec(f1,-4,0,40,1e-5);
  disp(sprintf('iter=%d, root=%d, residual=%d', n, x_bisec, fx));
  [n,x_bisec,fx]=bisec(f1,0,4,40,1e-5);
  disp(sprintf('iter=%d, root=%d, residual=%d', n, x_bisec, fx)); 
  
  
  f2=inline('exp(x)-2-cos(exp(x)-2)');
  df2=inline('exp(x)+exp(x)*sin(exp(x)-2)');
  [n,x_bisec,fx]=bisec(f2,-4,4,40,1e-5);
  disp(sprintf('iter=%d, root=%d, residual=%d', n, x_bisec, fx));

  
  f3=inline('exp(x)+2^(-x)+2*cos(x)-6');
  df3=inline('exp(x)+2^(-x)*log(2)-2*sin(x)');
  [n,x_bisec,fx]=bisec(f3,-4,0,40,1e-5);
  disp(sprintf('iter=%d, root=%d, residual=%d', n, x_bisec, fx));
  [n,x_bisec,fx]=bisec(f3,0,4,40,1e-5);
  disp(sprintf('iter=%d, root=%d, residual=%d', n, x_bisec, fx));

  fprintf('Ans. to Q(2)\n')
  disp('hit ENTER to proceed')
  pause
  
  %[n,x_newton,fx,X]=newton(f1,df1,-2,20,1e-5);
  %disp(sprintf('iter=%d, x0=%d, root=%d, residual=%d', n, -2, x_newton, fx));
  
  [n,x_newton,fx,X]=newton(f1,df1,2,20,1e-5);
  disp(sprintf('iter=%d, x0=%d, root=%d, residual=%d', n, 2, x_newton, fx));

  [n,x_newton,fx,X]=newton(f2,df2,1.5,20,1e-5);
  disp(sprintf('iter=%d, x0=%d, root=%d, residual=%d', n, 1.5, x_newton, fx));
  
  %[n,x_newton,fx,X]=newton(f3,df3,-3.5,20,1e-5);
  %disp(sprintf('iter=%d, x0=%d, root=%d, residual=%d', n, -3.5, x_newton, fx));
  
  [n,x_newton,fx,X]=newton(f3,df3,2,20,1e-5);
  disp(sprintf('iter=%d, x0=%d, root=%d, residual=%d', n, 2, x_newton, fx));
  
  fprintf('Ans. to Q(3)\n')
  disp('hit ENTER to proceed')
  pause

  n=50; tol=1e-10;
  x0=4;
  [n,xx,fx,Xn1]=newton(f1,df1,x0,n,tol);
  [n,x_bisec,fx,Xb1]=bisec(f1,0,4,50,tol);

  x0=3;
  [n,xx,fx,Xn2]=newton(f2,df2,x0,n,tol);
  [n,x_bisec,fx,Xb2]=bisec(f2,0,4,50,tol);

  x0=2;
  [n,xx,fx,Xn3]=newton(f3,df3,x0,n,tol);
  [n,x_bisec,fx,Xb3]=bisec(f3,0,4,50,tol);




  figure(1);clf
  semilogy(Xb1(:,1),Xb1(:,3),'o-')
  hold on
  semilogy(Xn1(:,1),Xn1(:,3),'d:')
  set(gca,'fontsize',16);
  xlabel('# of iterations');
  ylabel('|f(x^*)');
  title('x+3*cos(x)-exp(x)');
  
  axis tight; axis square
  legend('bisection','Newton')
  
  figure(2);clf
  semilogy(Xb2(:,1),Xb2(:,3),'o-')
  hold on
  semilogy(Xn2(:,1),Xn2(:,3),'d:')
  set(gca,'fontsize',16);
  xlabel('# of iterations');
  ylabel('|f(x^*)');
  axis tight; axis square
  legend('bisection','Newton')
  title('exp(x)-2-cos(exp(x)-2)');
  
  figure(3);clf
  semilogy(Xb3(:,1),Xb3(:,3),'o-')
  hold on
  semilogy(Xn3(:,1),Xn3(:,3),'d:')
  set(gca,'fontsize',16);
  xlabel('# of iterations');
  ylabel('|f(x^*)');
  axis tight; axis square
  legend('bisection','Newton')
  title('exp(x)+2^(-x)+2*cos(x)-6');
  
  return
  xmin=-1;
  xmax=1;
  
  for k=1:3
    %Tabulated data
    x=xmin:1/2^k:xmax;
    y=x.^2.*exp(-x./2);
    
    %Interpolating nodes
    xi=xmin:0.01:xmax;
    yi=xi.^2.*exp(-xi./2);
    
    [p,q,err]=Interpolate(x,y,xi,yi);
    figure(1); clf
    set(gca,'fontsize',16)
    plot(xi,p,'-',xi,q,'--',xi,yi,':')
    
    xlabel('x'); ylabel('p, q, y')
    legend('y','p','q')
    title(['h=' num2str(1/2^k)])
    
    figure(2); clf
    set(gca,'fontsize',16)
    plot(xi,err,'-');
    xlabel('x'); ylabel('error')
    title(['h=' num2str(1/2^k)])
    
    disp('hit ENTER to proceed')
    pause
  end

  fprintf('Ans. to Q(2)\n')
  disp('hit ENTER to proceed')
  figure(1);clf
  figure(2);clf
  pause
  xmin=-1;
  xmax=1;
  
  for k=1:3
    %Tabulated data
    x=xmin:1/2^k:xmax;
    y=sin(3.0*pi*x);
    
    %Interpolating nodes
    xi=xmin:0.01:xmax;
    yi=sin(3.0*pi*xi);
    
    [p,q,err]=Interpolate(x,y,xi,yi);
    figure(1); clf
    set(gca,'fontsize',16)
    plot(xi,p,'-',xi,yi,':')
    
    xlabel('x'); ylabel('p,  y')
    legend('p','y')
    title(['h=' num2str(1/2^k)])
    
    figure(2); clf
    set(gca,'fontsize',16)
    plot(xi,err,'-');
    xlabel('x'); ylabel('error')
    title(['h=' num2str(1/2^k)])
    
    disp('hit ENTER to proceed')
    pause
  end

  fprintf('Ans. to Q(3): y=sin(6pix)\n')
  disp('hit ENTER to proceed')
  figure(1);clf
  figure(2);clf
  pause
  xmin=-1;
  xmax=1;
  
  for k=1:3
    %Tabulated data
    x=xmin:1/2^k:xmax;
    y=sin(6.0*pi*x);
    
    %Interpolating nodes
    xi=xmin:0.01:xmax;
    yi=sin(6.0*pi*xi);
    
    [p,q,err]=Interpolate(x,y,xi,yi);
    figure(1); clf
    set(gca,'fontsize',16)
    plot(xi,p,'-',xi,yi,':')
    
    xlabel('x'); ylabel('p, y')
    legend('p','y')
    title(['h=' num2str(1/2^k)])
    
    figure(2); clf
    set(gca,'fontsize',16)
    plot(xi,err,'-');
    xlabel('x'); ylabel('error')
    title(['h=' num2str(1/2^k)])
    
    disp('hit ENTER to proceed')
    pause
  end
  return
  
function [yp,q,err]=Interpolate(x,y,xi,yi)
  
%Newton polynomial
  D=divdiff(x,y);
  p=newton_pol(xi,x,y);
  
  %Benchmark with matlab
  z=interp1(x,y,xi,'cubic');
  
  yp=p;
  q=z;
  err=abs(p-yi);
  
