% Input:
% A0 - initial conditions vector: [Kc vc y Db eta lambda] or [Kc vc Ks vs y Db eta lambda]
% flags - vector of settings: [realization matrix solver lf color]
%         realization: 0 (static); 1 (rigid); 2 (gauss); 3 (floppy);
%         matrix: 1 (spin symmetric); 2 (spin symmetric flow backwards); 3 (full); 4 (full backwards);
%         solver: 0 (explicit Runge-Kutta (4,5) formula, the Dormand-Prince pair)
%                 1 (variable order Adams-Bashforth-Moulton PECE solver)
%                 2 (modified Rosenbrock formula of order 2)
%                 3 (trapezoidal rule)
%                 4 (explicit Runge-Kutta (2,3) pair of Bogacki and Shampine)
%         lf: length of integration;
%         color: 0 (no plots); 1 (black); 2 (blue); 3 (red); 4 (green); 5 (magenta);
                 
function Af = flow(A0,flags,filename);

clrs = ['k' 'b' 'r' 'g' 'm'];
flowMatrix = matrixFileName(flags(1),flags(2));
lspan = 0:0.1:flags(4);
options = odeset('MaxStep',20000);%,'AbsTol',22,'RelTol',12);

switch flags(3)
    case 0
        [l,A]=ode45(flowMatrix, lspan, A0, options);
    case 1
        [l,A]=ode113(flowMatrix, lspan, A0);
    case 2
        [l,A]=ode23s(flowMatrix, lspan, A0);
    case 3
        [l,A]=ode23t(flowMatrix, lspan, A0, options);
    case 4
        [l,A]=ode23(flowMatrix, lspan, A0);
end

if strcmp(filename,'n') == 0
	dlmwrite(filename,[l A],' ');
end

if flags(5)~=0
    switch flags(1)
        case 0
            statPlot(l,A,clrs(flags(5)));
        case 1
            rigdPlot(l,A,clrs(flags(5)));
        case 2
            gausPlot(l,A,clrs(flags(5)));
        case 3
            flopPlot(l,A,clrs(flags(5)));
    end
end
Af = A(end,:);