% Setup: solve y' = f(t,y(t)), y(0) = y0, t <= T, with
g = @(t)(t.^4 - 6*t.^3 + 12*t.^2 - 14*t + 9)./(1+t).^2;
f = @(t,y)(y.^2 - g(t));
y0 = 2;
T = 2.2;

% Exact solution
yexact = @(t)(1-t).*(2-t)./(1+t);
hfine = 1e-3;
tfine = 0:hfine:T;
yfine = yexact(tfine);
figure(1); clf;
plot(tfine,yfine,'r-','linewidth',2)
axis([0 T -1 2]);
set(gca,'fontsize',16);

%% Numerical solution: Forward Euler (FE)

for h = [.2 .1 .05]
    t = 0:h:T;
    nsteps = length(t)-1;
    yFE = [y0; zeros(nsteps, 1)];    
    for n = 1:nsteps
        yFE(n+1) = yFE(n) + h*f(t(n),yFE(n));        
    end
    figure(1); hold on;
    plot(t,yFE,'b-x','linewidth',1);
    axis([0 T -1 2]);
end
xlabel('t'); ylabel('y'); title('Forward Euler')
legend('Exact','h = .2', 'h = .1', 'h = .05');

%% Numerical solution: Runge-Kutta 2 (RK)

figure(1); clf;
plot(tfine,yfine,'r-','linewidth',2)
set(gca,'fontsize',16);
for h = [.2 .1 .05]
    t = 0:h:T;
    nsteps = length(t)-1;
    yRK = [y0; zeros(nsteps, 1)];    
    for n = 1:nsteps
        yFE = yRK(n) + h*f(t(n),yRK(n));
        yRK(n+1) = yRK(n) + h/2*...
            ( f(t(n),yRK(n)) + f(t(n+1), yFE) );        
    end
    figure(1); hold on;
    plot(t,yRK,'b-x','linewidth',1);
    axis([0 T -1 2]);
end
xlabel('t'); ylabel('y'); title('Runge-Kutta 2')
legend('Exact','h = .2', 'h = .1', 'h = .05');

%% Numerical solution: Backward Euler (BE)

figure(1); clf;
plot(tfine,yfine,'r-','linewidth',2)
set(gca,'fontsize',16);
for h = [.1 .05]
    t = 0:h:T;
    nsteps = length(t)-1;
    yBE = [y0; zeros(nsteps, 1)];    
    for n = 1:nsteps
        F = @(x)(x - yBE(n) - h*f(t(n+1),x));        
        yBE(n+1) = fzero(F,yBE(n));
    end
    figure(1); hold on;
    plot(t,yBE,'b-x','linewidth',1);
    axis([0 T -1 2]);
end
xlabel('t'); ylabel('y'); title('Backward Euler')
legend('Exact','h = .1','h = .05');

%% Numerical solution: Midpoint method (MM)

figure(1); clf;
plot(tfine,yfine,'r-','linewidth',2)
set(gca,'fontsize',16);
for h = [.2 .1 .05]
    t = 0:h:T;
    nsteps = length(t)-1;
    yMM = [y0; zeros(nsteps, 1)];    
    for n = 1:nsteps
        F = @(x)(x - yMM(n) - h/2*...
            ( f(t(n),yMM(n)) + f(t(n+1),x)) );        
        yMM(n+1) = fzero(F,yMM(n));
    end
    figure(1); hold on;
    plot(t,yMM,'b-x','linewidth',1);
    axis([0 T -1 2]);
end
xlabel('t'); ylabel('y'); title('Midpoint method')
legend('Exact','h = .2', 'h = .1','h = .05');