N = 40;
t0 = 0;
tfin = 2;
h = (tfin - t0) / N;
tjs = (t0:h:tfin)';
y = [2/5];

% Define relevant functions
function out = f(t, yin)
  out = t.*yin.^2;
end

function out = fy(t, yin)
  out = 2*t.*yin;
end

function root = nrm(f, df, yLatest, tnext, h, tol, maxIt)
  z = yLatest;
  zold = z + 100;    % make z's initially quite far apart so loop will start
  k = 0;
  while ((norm(z - zold) >= tol) & (k < maxIt))
    k++;
    zold = z;
    z = zold - (yLatest - zold + h*f(tnext, zold)) / (h*df(tnext, zold) - 1);
  end
  root = z;
end

% solution process
for jj = 1:N
  ynew = nrm(@f, @fy, y(end), tjs(jj+1), h, 10^(-4), 10);
  y = [y; ynew];
end

ts = (0:.01:2)';
plot(ts, 2./(5 - ts.^2), 'b-')
hold on
plot(tjs, y, 'r*')
hold off