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 care about boundary 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);
                act_input = P.Input(t);
                for c=-1:1
                    for d=-1:1
                        if strcmp(P.Boundary, '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)=act_input(a+c,b+d);
                                stateregion(c+2,d+2)=x(a+c,b+d);
                            end
                        elseif strcmp(P.Boundary, '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)=act_input(inda,indb);
                            stateregion(c+2,d+2)=x(inda,indb);
                        elseif strcmp(P.Boundary, '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)=act_input(inda,indb);
                            stateregion(c+2,d+2)=x(inda,indb);
                        else
                            disp('Wrong boundary condition')
                        end
                    end
                end
            else
                inputregion=act_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

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

end
