function [ x, sol ] = Schnakenberg( n )
%SCHNAKENBERG Solve 2D Schnakenberg equation on a square
% Parameters and initial conditions taken from Ruuth (1995)
% (initial conditions slightly modified)
%
h=1/n;
lambda = 1000.;
du = 1.0;
dv = sqrt(10);
%
x = (0:n-1)/n; %linspace(0,1,n);
y = x;  % Tensor product grid
[A,B] = compactlaplacian(n);
% 
ex = ones(n)*diag(x);
ey = diag(y)*ones(n);
cxy = cos(2*pi*(ex+ey));
cxy2= sin(2*pi*(2*ex-ey));
u0 = 0.919145*ones(n)+0.0016*cxy;
v0 = 0.937903*ones(n)+0.0316*cxy2;
for j=1:8,
    c = cos(2*pi*j*ex);
    u0 = u0 + 0.01*c;
    v0 = v0 + 0.01*c';
end
u0 = u0';
v0 = v0';
inif = [u0(:);v0(:)]; % Row order
% Compact fourth-order space derivatives
% Mehrstellen verfahren (Collatz, 1950)

% Mass matrix [[A,0],[0,A]]
 [i,j,s] = find(A);
 M = sparse([i,i+n^2],[j,j+n^2],[s,s],2*n^2,2*n^2);
       
    function yp = Sc(t,y)
        u = y(1:n^2,:);
        v = y(n^2+1:end,:);
        yp = zeros(size(y));
        usv = u.^2.*v;
        yp(1:n^2,:) = du^2*B*u/h^2+A*lambda*(0.126779 - u + usv);
        yp(n^2+1:end,:)= dv^2*B*v/h^2+A*lambda*(0.792366 - usv);
    end

opts = odeset('reltol',1.0e-7,'Mass',M,'Vectorized','on');
ts = 2.0;
sol = ode15s(@Sc,[0,ts], inif, opts );

end

