function [I,S] = tsw(I0,S0,N,sg)
% Taylor series coeffs of solution of S' = -S^3*I, I'= -S*I^2
% S(tn)=S0, I(tn)=I0; sg is the direction of integration.
% if tn=0, I0=S0=1, then I=exp(-W(t)) and S = 1/(1+W(t)).
    I = zeros(1,N+1);
    S = I; C = I; D = I; E = I; F = I;
    I(1) = I0; S(1) = S0;
    for k=1:N,
        C(k) = S(1:k)*I(k:-1:1).';
        D(k) = C(1:k)*I(k:-1:1).';
        E(k) = S(1:k)*S(k:-1:1).';
        F(k) = C(1:k)*E(k:-1:1).';
        S(k+1) = -sg*F(k)/k;
        I(k+1) = -sg*D(k)/k;
    end  
end