First you can transform the matrix to a graph, where every maximal area containing the same integer is transformed to a vertex with a weight that is the number of cells in the area, and there is an edge between two vertices iff the two corresponding areas are adjacent, that is, there is a cell in one area that is adjacent to at least one cell in the other area. Therefore, the matrix in your example is transformed to

    5 - 3 - 2 - 5
    |   | / |   | \
    2 - 1   |   |   3
    | / | \ |   |   |
    6   |   5 - 2 - 5
    |   | / | \ | / |
    5 - 2   |   3   |
    |   |   |   |   |
    4 - 6 - 8 - 9 - 6

Now the problem turns out to be finding the connected subgraph with maximum sum of weights such that it contains exactly two integers. This transformation can be done in $O(mn)$ time.

Next, run the following algorithm:

    maxS := empty set
    While True:
        Search the next edge (u,v) that is not used
        If not found:
            return maxS
        Let x,y be the integers contained in u,v respectively
        (#) Do BFS from u where only edges between vertices containing x,y respectively are considered
        Mark all edges found during the BFS as used
        S := the set of vertices found during the BFS
        If the sum of weights in S > the sum of weights in maxS
            maxS := S

Note the (#) line, this algorithm requires us to, given a vertex containing integer `x`, efficiently find all its neighbors (as well as edges) containing `y`. You can sort all the neighbors for all vertices in advance, which takes $O(\sum_v n_v\log n_v)=O(mn\log(mn))$ time ($n_v$ is the number of neighbors of vertex $v$). Now compared to normal BFS, it takes extra $O(\log n_v)$ time for each vertex $v$.

Basically, each BFS induces a maximal connected subgraph containing exactly two integers. This algorithm checks all such subgraphs and choose the optimal one. Note each edge belongs to only one such subgraph and each vertex $v$ belongs to $t_v$ such subgraphs ($t_v$ is the number of distinct integers among all neighbors of $v$), the time complexity of this algorithm is 

$$O\left(mn+\sum_vt_v\log n_v\right)=O\left(mn\log(mn)\right).$$