function s4(alpha,beta,gamma1,gamma2,T)

% ======================================================================
% Input:  T: ending time of simulations
% ======================================================================


  if ( nargin < 1 )
    alpha = 0.044;
    beta = 0.065;
    gamma1 = 2;
    gamma2 = 2;
    T = 0.050;
  end

% Define model parameters: MAKE SURE TO SET THE PARAMTERS CORRECTLY!!
D1 = 0.01;           %1;
D2 = 1;              %8.66;

% Define numerical parameters
k = 0.005;
t = 0:k:T;
%ckeep = 100000;             % MAKE SURE TO SET CKEEP CORRECTLLY.

L = 8;
Nx = 1024;                 % N intervals in space;
Ny = 1024;
Nx = 256;                 % N intervals in space;
Ny = 256;
x = (0:L/Nx:L)';           % Grid points in physical space;
y = x;
[xx,yy]=meshgrid(x,y);
u0= beta;
v0= (alpha+beta)/(2*beta^2);

%_________________________________________

Hatim=fopen('x.dat', 'w'); fopen(Hatim);
xcon=sprintf('%f ',x');fprintf(Hatim, '%s\n',xcon);
fclose(Hatim);
Hatim=fopen('y.dat', 'w'); fopen(Hatim);
ycon=sprintf('%f ',y');fprintf(Hatim, '%s\n',ycon);
fclose(Hatim);
% =====================================================================
[u, v] = initialize(xx, yy, Nx, Ny,u0,v0);

%tkeep = zeros(1, 1+(length(t)-1)/ckeep);
%ukeep = zeros(Nx, Ny, 1+(length(t)-1)/ckeep);
%vkeep = zeros(Nx, Ny, 1+(length(t)-1)/ckeep);
%_________________________________
Hatim1=fopen('t.dat','w'); fopen(Hatim1);
Hatim2=fopen('u.dat','w'); fopen(Hatim2);
Hatim3=fopen('v.dat','w'); fopen(Hatim3);
fprintf(Hatim1, '%4.8f\n',t(1));
for n=1:Nx
ucon=sprintf('%f ',u(n,:)); fprintf(Hatim2, '%s\n',ucon);
vcon=sprintf('%f ',v(n,:)); fprintf(Hatim3, '%s\n',vcon);
end
% tkeep(1)   = 0;
% ukeep(:,:,1) = u;
% vkeep(:,:,1) = v;

sx=((2*pi/L)*[0:Nx/2 -Nx/2+1:-1]');
sy=((2*pi/L)*[0:Ny/2 -Ny/2+1:-1]');
[sxx,syy] = meshgrid(sx, sy);
s1 = (sxx.^2+syy.^2).^(gamma1/2);
s2 = (sxx.^2+syy.^2).^(gamma2/2);

f = @(t,uv)schnak_deriv(t,uv,D1,D2,s1,s2,alpha,beta,Nx,Ny);
tspan = [ 0.0, T ];
u = reshape ( u, Nx * Ny, 1 );
v = reshape ( v, Nx * Ny, 1 );
uv0 = [ u; v ];
size ( uv0 )
dt = 0.005;
tol = 1.0E-2;
%[t, uv, e ] = rk12_adapt ( f, tspan, uv0, dt, tol );
options = odeset ( 'reltol', 1.0e-2, 'abstol', 1.0e-2 );
[ t, uv ] = ode45 ( f, tspan, uv0, options );

length ( t )

% =====================================================================
%Initial Condition 
function [ut, vt] = rhsfun(u, v, s1, s2, alpha, beta)

U = fft2(u);
V = fft2(v);

Lu = ifft2(U*s1);
Lv = ifft2(V*s2);

ut = Lu+funf(u,v,alpha,beta);
vt = Lv+fung(u,v,alpha,beta);

function uvt = schnak_deriv (t, uv, D1, D2, s1, s2, alpha, beta,Nx,Ny)

u = uv(1:Nx*Ny);
u = reshape ( u, Nx, Ny );
v = uv(Nx*Ny+1:2*Nx*Ny);
v = reshape ( v, Nx, Ny );

U = fft2(u);
V = fft2(v);

Lu = real ( ifft2(-D1*U*s1) );
Lv = real ( ifft2(-D2*V*s2) );

ut = Lu+funf(u,v,alpha,beta);
ut = reshape ( ut, Nx*Ny, 1 );
vt = Lv+fung(u,v,alpha,beta);
vt = reshape ( vt, Nx*Ny, 1 );
uvt = [ ut; vt];

fprintf ( 1, 't = %g\n', t );

% =====================================================================
%Initial Condition 
function [u, v] = initialize(xx, yy, Nx, Ny,u0,v0)
u = u0*ones(Nx,Ny);
v = v0*ones(Nx,Ny);

m = (Nx-20)/2;
n = (Ny-20)/2;
u(m:m+20, n:n+20) = 0.5;
v(m:m+20, n:n+20) = 0.25;

% =====================================================================
function [output] = funf(u,v,alpha,beta)
output =((beta-alpha)/2)+u.^2.*v;

function [output] = fung(u,v,alpha,beta)
output = ((beta+alpha)/2)-u.^2.*v;
