function [xx,c,t_out,a1,a2,a3,a4,a5]=ch_oned_solve(t_final)

% *************************************************************************
% One-diemnsional pseduospectral code to solve the Cahn-Hilliard equation
%
% \partial_t C = D\partial_{xx}(C^3-C-\gamma \partial_{xx}C), t>0,
% x\in(0,L),
%
% with periodic boundary conditons C(x+L,t)=c(x,L), t>0, and initial
% conditions C(x,t=0)=C_{init}(x), x\in [0,L].
%
% D is hardcoded to be one.
%
% *************************************************************************
% Numerical parameters

[L,gamma,Nx,dt,n_max,my_phases,ampl_init,ampl_pert]=fix_all_parameters();

% Fundamental wavenumber:
k0=2*pi/L;

% Grid spacing:
dx=L/Nx;

% Spatial grid:
xx=(0:(Nx-1))*dx;

% Wavenumber grid:
kx=(2*pi/L)*(-(Nx/2):1:((Nx/2)-1));

% Final timestep:
t_ctr_final=floor(t_final/dt);

% Vector of timesteps:
t_out=(1:t_ctr_final)*dt;

% *************************************************************************
% Initial conditions - made of of fundamental plus overtones.  The
% overtones are stored in the array pert.  Pert contains n_max overtones.

pert=0*xx;

for n=2:n_max
    pert=pert+cos(n*k0*xx+2*pi*my_phases(n));
end

pert=pert/(n_max-2+1);

c_init=ampl_init*( cos(k0*xx) + ampl_pert*pert );
c=c_init;

% *************************************************************************
% Arrays to store Fourier amplitudes:

a1=0*(1:t_ctr_final);
a2=0*(1:t_ctr_final);
a3=0*(1:t_ctr_final);
a4=0*(1:t_ctr_final);
a5=0*(1:t_ctr_final);

% *************************************************************************
% Step into time loop

for t_ctr=1:t_ctr_final
    
    % PDE code - pseudospectral method.
    
    mu=c.^3-c;
    mu_hat=fftshift(fft(mu));
    c_hat=fftshift(fft(c));
        
    RHS=c_hat+dt*(-kx.*kx.*mu_hat);
    c_hat=RHS./(1+gamma*dt*kx.*kx.*kx.*kx);
    
    c=ifft(ifftshift(c_hat));
    
    % (Nx/2)+1 is the k=0 mode.
    % (Nx/2)+2 is the k=(2*pi/L) mode.
    % Thereafter it is multiples of (2*pi/L).
    
    % Here we store the amplitudes of the Fourier modes as generated by the
    % PDE pseudospectral method.
    
    a1(t_ctr)=c_hat((Nx/2)+2);
    a2(t_ctr)=c_hat((Nx/2)+3);
    a3(t_ctr)=c_hat((Nx/2)+4);
    a4(t_ctr)=c_hat((Nx/2)+5);
    a5(t_ctr)=c_hat((Nx/2)+6);
    
end

% *************************************************************************


end