function [t, y] = solutions_Final_p1c(N)

    t = linspace(0, 1, N);
    y = zeros(size(t));
    
    dt = t(2) - t(1);
    
    for i=2:length(t)
        k1 = solutions_Final_p1a(t(i-1), y(i-1));
        k2 = solutions_Final_p1a(t(i), y(i-1) + dt*k1);
        y(i) = y(i-1) + 0.5*dt*(k1 + k2);
    end
    

end

