function dx = cell_equation_boundary(t,X,P)

%reshape the 1xN input for the size of the image
x=reshape(X,P.SizeX,P.SizeY);

dx=zeros(P.SizeX,P.SizeY);

%we do not care about boundar conditions now
for a=1:P.SizeX
    for b=1:P.SizeY
	if (a==1) || (b==1) || (a==P.SizeX) || (b==P.SizeY)
		inputregion=zeros(3,3);
		stateregion=zeros(3,3);
		for c=-1:1
			for d=-1:1
			    if P.Bounadry=='Constant'
				if a+c<1 | b+d<1 | a+c>P.SizeX | b+d>P.SizeY
					inputregion(c+2,d+2)=P.BounValue;
					stateregion(c+2,d+2)=P.BounValue;
				else
					inputregion(c+2,d+2)=P.Input(a+c,b+d);
					stateregion(c+2,d+2)=x(a+c,b+d);
				end
			    elseif P.Bounadry=='ZeroFlux'
			    	inda=a+c;
				if a+c<1
					inda=1;
				elseif a+c>P.SizeX
					inda=P.SizeX;   	
				end
			    	indb=b+d;
				if b+d<1
					indb=1;
				elseif b+d>P.SizeY
					indb=P.SizeY;   	
				end
				inputregion(c+2,d+2)=P.Input(inda,indb);
				stateregion(c+2,d+2)=x(inda,indb);
			    elseif P.Bounadry=='Periodic'
			    	inda=a+c;
				if a+c<1
					inda=P.SizeX;
				elseif a+c>P.SizeX
					inda=1;   	
				end
			    	indb=b+d;
				if b+d<1
					indb=P.SizeY;
				elseif b+d>P.SizeY
					indb=1;   	
				end
				inputregion(c+2,d+2)=P.Input(inda,indb);
				stateregion(c+2,d+2)=x(inda,indb);
			    else
				disp('Wrong boundary condition')
			    end
			end
		end
	else
		inputregion=P.Input(a-1:a+1,b-1:b+1);	    
		stateregion=x(a-1:a+1,b-1:b+1);
		
	end
	u=reshape(inputregion,9,1);
	stateregion=reshape(stateregion,9,1);
        y=nonlinearity1(stateregion);

        dx(a,b)=-x(a,b) + P.A*y  + P.B*u + P.Z;
    end
end

dx(1,1)=dx(a,b)+4.04*sin(pi*t/2)

%reshape back fro Nx1
dx=reshape(dx,P.SizeX*P.SizeY,1);