function [t, y] = solutions_hw7_p4(f, tspan, y0, order, N)
    switch order
        case 1
            fcn = @FEuler;
        case 2
            fcn = @IEuler;
        case 4
            fcn = @RK4;
        otherwise
            error('Unknown order');
    end
    
    t = linspace(tspan(1), tspan(2), N);
    dt = t(2) - t(1);
    y = zeros(size(t));
    y(1) = y0;    
    for i=2:N
       y(i) = fcn(f, t(i-1), y(i-1), dt);
    end
end


function [y] = FEuler(f, t, y, dt)
    k = f(t, y);
    y = y + dt*k;
end


function [y] = IEuler(f, t, y, dt)
    k1 = f(t, y);
    k2 = f(t + dt, y + dt*k1);
    y = y + 0.5*dt*(k1 + k2);
end

function [y] = RK4(f, t, y, dt)
    k1 = f(t, y);
    k2 = f(t + 0.5*dt, y + 0.5*dt*k1);
    k3 = f(t + 0.5*dt, y + 0.5*dt*k2);
    k4 = f(t + dt, y + dt*k3);
    
    y = y + dt*(k1 + 2.0*k2 + 2.0*k3 + k4)/6.0;
end