clear all
close all
a = [0,0;1/2,1/2];
c(2) = 1;
b = [1/2,1/2];
y0 = 1;
tsrif = 100;
ts = tsrif;
y = y0;
k = 1/ts;
t = 0;
g = @(t,y) t*sin(y);
dg = @(t,y) t*cos(y);
tol = 1e-5;
for n = 1:ts
  f1 = 2*y-g(t,y);
  F = @(f2) f2-2*(y+k*a(2,1)*f1+k*a(2,2)*f2)-...
      g(t+c(2)*k,y+k*a(2,1)*f1+k*a(2,2)*f2);
  J = @(f2) 1-2*k*a(2,2)-dg(t+c(2)*k,y+k*a(2,1)*f1+k*a(2,2)*f2)*k*a(2,2);
  f2 = 2*y-g(t,y);
  errest = -J(f2)\F(f2);
  while (abs(errest) > tol)
    f2 = f2+errest;
    errest = -J(f2)\F(f2);
  end
  y = y+k*b(1)*f1+k*b(2)*f2;
  t = t+k;
end
yrif = y;
tsrange = 10:10:50;
counter = 0;
for ts = tsrange
  counter = counter+1;
  y = y0;
  k = 1/ts;
  t = 0;
  g = @(t,y) sin(y)*t;
  dg = @(t,y) cos(y)*t;
  tol = 1e-5;
  for n = 1:ts
    f1 = 2*y-g(t,y);
    F = @(f2) f2-2*(y+k*a(2,1)*f1+k*a(2,2)*f2)-...
	g(t+c(2)*k,y+k*a(2,1)*f1+k*a(2,2)*f2);
    J = @(f2) 1-2*k*a(2,2)-dg(t+c(2)*k,y+k*a(2,1)*f1+k*a(2,2)*f2)*k*a(2,2);
    f2 = 2*y-g(t,y);
    errest = -J(f2)\F(f2);
    while (abs(errest) > tol)
      f2 = f2+errest;
      errest = -J(f2)\F(f2);
    end
    y = y+k*b(1)*f1+k*b(2)*f2;
    t = t+k;
  end
  errore(counter) = abs(yrif-y);
end
loglog(tsrange,errore,tsrange,tsrange.^(-2))