prob4ptC.m

function x = prob4ptC(prob4_ode, dt, x0,tlimits)
    t = [tlimits(1):dt:tlimits(2)];
    x = zeros(length(t),length(x0));
    x(1,:) = x0;
    for i = 2:length(t)
        dx = prob4_ode(t(i-1),x(i-1,:));
        x(i,:) = x(i-1,:)+dt*dx;
        dxc = (dx + prob4_ode(t(i),x(i,:)))/2;
        x(i,:) = x(i-1,:)+dt*dxc;
        n = 1;
        while (1) 
            n = n+1;
            yold = x(i,:);
            dxc = (dx+prob4_ode(t(i),x(i,:)))/2;
            x(i,:) = x(i-1,:)+dt*dxc;
            ea = abs(x(i,:)- yold)./x(i,:)*100; 
            if max(ea) < 0.0001 || n>100; break, end
        end
    end
end