function [unew] = stepRD(u,dt,Diff)
%
%   This routine updates the PDE variable u for 
%      u_t = u_{xx}+u-u^3, 
%       u(0,x)=0.5,          u_x(t,0)=u_x(t,1)=0
%
%   Inputs:    u - values at the current time
%              dt- the size of the timestep
%              Diff- The Diffusion coefficient (a scalar constant)
%
%   Outputs:   unew -values at the new time
n=length(u);
%%%%%%%%%% Initialize
% set up h   
h=(1-0)/(n-1);

% One on diagonal except 1st and last row
rows=2:n-1;
ID=sparse(rows,rows,1,n,n);
% Ones above the diagonal
UP=sparse(rows,rows+1,1,n,n);
% Ones below the diagonal
DN=sparse(rows,rows-1,1,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);

L=(UP-2*ID+DN)/h^2;
%Now form left hand side matrix
LHS=ID-dt.*Diff.*L/2 + bc;
RHS=(ID+dt.*Diff.*L/2)*u + dt.*(u-u.^3);
%
RHS(1)=0;
RHS(n)=0;

% equation is LHS*u = RHS, so invert LHS:
unew = LHS\RHS;
