function [u]=SOLVEGL(N,X,n,T,method,plotgraph);
%
% Solves the discretized version of the complex Ginzburg-Landau 
% equation: 
%
% u_t = alpha u_xx + epsilon*u -beta*|u|^2*u    
%
% with:  1. x in [-X,X] 
%        2. t in [0,T]
%        3. periodic boundary conditions 
%        4. u0(x)=0.8*(1/cosh(x-X/10)^2+1/cosh(x+X/10)^2)
%
% N : number of space points
% n : number of time steps
%
% If plotgraph=1 the solution is plot 
%
% Parameters of Ginzburg-Landau
%
epsilon=1;
c1=1;
c3=-2;
alpha=1+i*c1;
beta=1-i*c3;
%
% A is the Discritized Laplacian
%
A=zeros(N,N);
%
% First and last lines of A
%
A(1,1)=-2;
A(1,2)=1;
A(1,N)=1;
%
A(N,1)=1;
A(N,N-1)=1;
A(N,N)=-2;
%
% Other lines
%
for p=2:N-1,
    A(p,p-1)=1;
    A(p,p)=-2;
    A(p,p+1)=1;
end;
%
Dx=2*X/(N+1);
A=1/Dx^2*A;
%
%
% Coefficients of the splitting method as in 
% exp(a_1 hA) exp(b_1 hB)...exp(a_s hA) exp(b_s hB) 
%
if method==1,
    %
    % Strang-Splitting
    %
    s=2;
    a=[1;0];
    b=[1/2;1/2];
    %
elseif method==2,
    %
    % Method of CCDV10'
    %
    s=5;
    a=[1/4;1/4;1/4;1/4;0];
    b=[1/10-i/30;4/15+2*i/15;4/15-i/5;4/15+2*i/15;1/10-i/30];
    %
elseif method==3,
    %
    % Method of Blanes, Casas, Chartier and Murua 
    %
    s=17;
    a=1/16*[1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;1;0];
    b=[0.028920177910074098791 - 0.005936580835725746103*i;
       0.056654351383649876160 + 0.020841963949772627119*i;
       0.067258385822722143569 - 0.039386393748812362460*i;
       0.070333980553260772061 + 0.058952097930307840316*i;
       0.077095100838099173580 - 0.038247636602014810025*i;
       0.042022140317231098258 - 0.033116379859951038579*i;
       0.050147397749937784280 + 0.061283684958324249562*i;
       0.047750191909146143447 - 0.032332468814362628262*i;
       0.119636547031757819706 + 0.015883426044923736862*i;
       0.047750191909146143447 - 0.032332468814362628262*i;
       0.050147397749937784280 + 0.061283684958324249562*i;
       0.042022140317231098258 - 0.033116379859951038579*i;
       0.077095100838099173580 - 0.038247636602014810025*i;
       0.070333980553260772061 + 0.058952097930307840316*i;
       0.067258385822722143569 - 0.039386393748812362460*i;
       0.056654351383649876160 + 0.020841963949772627119*i;
       0.028920177910074098791 - 0.005936580835725746103*i];
end;
%
h=T/n;
if plotgraph==1, 
    us=zeros(N,n);
end;
%
% Initial values
%
u=zeros(N,1);
x=zeros(N,1);
for p=1:N,
    x(p)=-X+p*Dx;
end;
u=0.8*sech((x+10).^2) + 0.8*sech((x-10).^2);
v=real(u);
w=imag(u);
%
time=zeros(n,1);
%
% In what follows, the ai's are assumed to be real!
%
EA(:,:,1)=exp(epsilon*h*a(1))*expm(h*alpha*a(1)*A);
EAb(:,:,1)=conj(EA(:,:,1));
for k=2:s,
    if a(k)==a(k-1),
        EA(:,:,k)=EA(:,:,k-1);
        EAb(:,:,k)=conj(EA(:,:,k));
    else
        EA(:,:,k)=exp(epsilon*h*a(k))*expm(h*alpha*a(k)*A);
        EAb(:,:,k)=conj(EA(:,:,k));
    end;
end;
%
% Main loop
%
for r=1:n,
    time(r)=r*h;
    vt=(-i/2*v+1/2*w);
    wt=(1/2*v-i/2*w);
    for k=1:s,
        %
        % One step of exp(b1(k)*h*B) corresponding to 
        % u'=-beta*|u|^2*u
        %
        Mt=4*i.*vt.*wt;
        L=log(1+2*h*b(k)*Mt);
        vt=exp(-beta/2*L).*vt;
        wt=exp(-conj(beta)/2*L).*wt;
        %
        % One step of exp(a1(k)*h*A) corresponding to 
        % u'=alpha * Delta + eps*u
        %
        vt=EA(:,:,k)*vt;
        wt=EAb(:,:,k)*wt;
    end;
    v=real(i*vt+wt);
    w=real(vt+i*wt);
    u=v+i*w;
    if plotgraph==1, 
        us(:,r)=u;
    end;
end;
if plotgraph==1, 
    figure(1);
    clf;
    surf(time,x,real(us),'EdgeColor','none');
    title('Solution of the Ginzburg-Landau equation (real part)');
    xlabel('Time');
    ylabel('x');
    zlabel('u(t,x)');
    view(90,90),shading interp
    colormap jet
%
    figure(2);
    clf;
    surf(time,x,imag(us),'EdgeColor','none');
    title('Solution of the Ginzburg-Landau equation (imaginary part)');
    xlabel('Time');
    ylabel('x');
    zlabel('u(t,x)');
    view(90,90),shading interp
    colormap jet
%
    figure(3);
    clf;
    surfc(time,x,abs(us),'EdgeColor','none');
    view(90,90),shading interp
    colormap jet;
    title('Solution of the Ginzburg-Landau equation (amplitude)');
    xlabel('Time');
    ylabel('x');
    zlabel('u(t,x)');
end;