% trapezoidal rule time stepping example 12.3 from Suli/Mayers
% Note that in the book the sign of g is incorrect
g = @(x) (x.^4-6*x.^3+12*x.^2-14*x+9)/((1+x).^2);
f = @(x,y) y.^2-g(x);
exact_y= @(x) (1-x).*(2-x)./(1+x);
x0 = 0;
y0 = 2;
% number of steps
N = 10;
% final time
T = 2;
% time step size
h = T/N;

method = 3; % 1 = modEuler; 2 = impEuler; 3 = RK4

X = zeros(1,N+1);
xk = x0;
Y = zeros(1,N+1);
yk = y0;
Y(1) = y0;
for k=1:N
   xk = x0 + (k-1)*h;
   if (method == 1)
       yk = yk + h*f(xk+h/2,yk+h/2*f(xk,yk));
   elseif (method == 2)
       k1 = f(xk,yk);
       yk = yk + h/2*(k1+f(xk+h,yk+h*k1));
   elseif (method == 3)
       % classical 4th order RK
       % compute stages--note that each requires evaluation of f
       k1 = f(xk,yk);
       k2 = f(xk+h/2,yk+k1*h/2);
       k3 = f(xk+h/2,yk+k2*h/2);
       k4 = f(xk+h,yk+h*k3);
       yk = yk + h/6*(k1+2*k2+2*k3+k4);
   else
       error('Must choose a method!')
   end
   % keep record for plotting
   X(k) = xk;
   Y(k+1) = yk;
end

X(N+1) = x0 + N*h;
plot(X,Y,'ro','linewidth', 2); hold on
Xfine = linspace(0,T,200);
plot(Xfine,exact_y(Xfine),'b','linewidth', 2)