%
%  script for hmwk #0:  convergence study
%
%  save as: "hw00c.m"
%
%- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

%  define (empty) data vectors

Nv = [];  Yv = [];  Ev = [];

%  set N loop (25 runs on my unix machine)

%nnmax = 25
nnmax = 16

for nn = 6:1:nnmax
  N = 2^nn;  disp(['the value of N = ' num2str(N)])
  h = pi/N;

%  execute recursion loop (no storage!)

  yold = 0;
  for j=2:1:N+1
    ynew = yold + h*(1 + yold^2)/4;
    yold = ynew;
  end

%  show the N+1th value (simple & fancy output)

  format long; disp(' ')
  disp(['the N+1th value of y = ' num2str(ynew,14)])
  disp(' '); format short

%  append new data

  Nv = [Nv N];
  Yv = [Yv ynew];
  Ev = [Ev abs(1 - ynew)];
end

clf 
loglog(Nv(2:end),Ev(2:end),'x')
hold on 
loglog(Nv(2:end),2*Ev(2:end)-Ev(1:end-1),'o')
title('\bf comparison of two limits')
xlabel('\bf N \rm (log-scale)')
ylabel('\bf error(N) \rm (log-scale)')