function gFlow(ini0,filename,clr,solver)

linc = 0.1;
lf = 8;
lspan = 0:linc:lf;
a = 0.1;
N = 4;
zN = 15;
L = (N-1)*a;
u = 0:a:L;
k = (1/N).*(0:N-1)';
disAmp = ini0(6);
%ini = vertcat(ini0(1:3),ini0(5:end),[zeros(zN,1); disAmp*ones(N-zN,1)]);
ini = vertcat(ini0(1:5),ini0(7:end),disAmp*ones(N,1));
Y = gaussEvalFlow(ini,lspan,solver);
SY = size(Y);

l = Y(:,1); Kc = Y(:,2); vc = Y(:,3); Ks = Y(:,4); vs = Y(:,5);
y = Y(:,6); eta = Y(:,7); lambda = Y(:,8);
Dk = Y(:,9:N+8); Du = Y(:,N+9:2*N+8);
[l(end) Kc(end) Dk(end)]

if strcmp(filename,'nf') == 0
    txt_a = strcat('%a=',num2str(a));
    txt_N = strcat('%N=',num2str(N));
    txt_lf = strcat('%lf=',num2str(lf));
    txt_linc = strcat('%linc=',num2str(linc));
    txt_sol = strcat('%solver=',num2str(solver));
    txt_ini = strcat('%ini=',num2str(ini0));
    header = strvcat(txt_a,txt_N,txt_lf,txt_linc,txt_sol,txt_ini);
    %dlmwrite(filename,header,' ');
	%dlmwrite(filename,Y,'-append','roffset',1,'delimiter',' ');
    dlmwrite(strcat('hdr',filename),header,' ');
    dlmwrite(filename,Y*1e3,' ');
end

if strcmp(clr,'np') == 0
    D1 = Du(1,:);
    D2 = Du(2,:);
    D3 = Du(3,:);
    Dk1 = Dk(1,:);
    Dk2 = Dk(2,:);
    Dk3 = Dk(3,:);

%     figure(3);
%     subplot(3,1,1)
%     plot(l,Dk(:,1),'b',l,Dk(:,round(0.5*N)),'g',l,Dk(:,N),'k');
%     grid on;
%     legend('1st component of D_{b}(k)',strcat(num2str(round(0.5*N)),'th component of D_{b}(k)'),strcat(num2str(N),'th component of D_{b}(k)'));
%     
%     subplot(4,1,2)
%     hold on;
%     plot(u,D1,'b');
%     plot(u,D5,'g');
%     plot(u,D10,'r');
%     title('D_{b}(u)');
%     xlabel('u');
%     ylabel('D');
%     grid on;
%     
%     subplot(4,1,3)
%     esthecis2(u,D1,'b');
%     esthecis2(u,D5,'g');
%     esthecis2(u,D10,'r');
%     
%     subplot(4,1,4)
%     hold on;
%     plot(k,Dk1,'b');
%     plot(k,Dk10,'g');
%     plot(k,Dk20,'r');
%     title('D_{b}(k)');
%     xlabel('k');
%     ylabel('D');
%     grid on;
    
    figure(1);
%     subplot(2,1,2)
    hold on;
    plot(l,lambda,clr);
    title('\lambda (l)');
    ylabel('\lambda');
    xlabel('l');
    grid on;

%     subplot(2,1,1)
%     hold on;
%     plot(l,eta,clr);
%     title('\eta (l)');
%     ylabel('eta');
%     xlabel('l');
%     grid on;
    
    figure(2)
    subplot(2,2,1)
    hold on;
    plot(l,Dk(:,1),clr);
    ylabel('D');
    xlabel('l');
    grid on;
    
    subplot(2,2,2)
    hold on;
    plot(l,Kc,clr);
    ylabel('Kc');
    xlabel('l');
    grid on;
    
    subplot(2,2,3)
    hold on;
    plot(l,y,clr);
    ylabel('y');
    xlabel('l');
    grid on;
    
    subplot(2,2,4)
    hold on;
    plot(l,Ks,clr);
    plot(l,kysu2(y),strcat('-.*',clr));
    ylabel('Ks');
    xlabel('l');
    grid on;
    
    figure(3);
    subplot(2,1,1)
    hold on;
    plot(Kc,Du(:,1),clr);
    title('K_{c}-D_{b}(u=0) phase space');
    xlabel('K_{c}');
    ylabel('D_{b}');
    grid on;
    
    subplot(2,1,2)
    hold on;
    plot(Kc,Dk(:,1),clr);
    title('K_{c}-D_{b}(k=0) phase space');
    xlabel('K_{c}');
    ylabel('D_{b}');
    grid on;
end
