function [xnew,ynew] = stepbruss(x,y,n,dt,A,B,D1,D2)
%
%   This routine updates the Brusselator variables
%   x and y and updates them one time step of size dt
%
%   Inputs:    x -  the first reactact concentration
%              y -  the second reactant concentration
%              n - the number of points in x and y
%              dt- the size of the timestep
%            A,B,D1,D2 - parameters in the model
%
%   Model:     X_t = D1 X_ss + A - (B+1)X + X^2*Y
%              Y_t = D2 Y_ss + B*X - X^2*Y
%        boundary conditions: X_s(0,t)=X_s(1,t)=Y_s(0,t)=Y_s(1,t)=0
%
%   Outputs:   xnew - the new first reactant concentration
%              ynew - the new second reactant concentration

%%%%%%%%%% Initialize
% set up h   *** note: n must be chosen before this is run ***
h=(1-0)/(n-1);

% One on diagonal except 1st and last row
ID=sparse(2:n-1,2:n-1,ones(1,n-2),n,n);
% Ones above the diagonal
UP=sparse(2:n-1,3:n,ones(1,n-2),n,n);
% Ones below the diagonal
DN=sparse(2:n-1,1:n-2,ones(1,n-2),n,n);

% first and last line for boundary conditions
%  These values are the coefficients of a three point derivative on the boundary
values=[-3/(2*h) 2/h -1/(2*h) -3/(2*h) 2/h -1/(2*h)];
bc=sparse([1 1 1 n n n],[1 2 3 n n-1 n-2],values,n,n);

%Now form left hand side matrix
LHSx=ID-dt.*D1.*(UP-2*ID+DN)./(2.*h^2) + bc;
LHSy=ID-dt.*D2.*(UP-2*ID+DN)./(2.*h^2) + bc;
RHSx=(ID+dt.*D1.*(UP-2*ID+DN)./(2.*h^2))*x + dt.*(A-(B+1).*x+x.^2.*y) ;
RHSy=(ID+dt.*D2.*(UP-2*ID+DN)./(2.*h^2))*y + dt.*(B.*x-x.^2.*y) ;
%
RHSx(1)=0;
RHSx(n)=0;
RHSy(1)=0;
RHSy(n)=0;

% the equation is LHS*x = RHS, so invert LHS:
xnew = LHSx\RHSx;
ynew = LHSy\RHSy;
