I don't know if the learning accuracy will be correct, but I will show the code that works without error.
The main questions in the questioner's code are as follows. I fixed that point.

For some reason, the pre-processing is limited to only two images.

In the model creation part, the input dimension of the site referenced by the questioner is completely different. It does not apply to the situation of the questioner. By the way, since the questioner is premised on images, I have included a CNN layer.

As I looked through the code, I referred to an example of a "merge model" in Keras. According to it, I think that the basic pattern is to merge completely different learning results such as images and text. I think that merging equivalent images and images like the questioner is essentially the same as synthesizing images into one and learning with a normal model.

from keras import backend as K
from keras.datasets import mnist
from keras.layers import Activation, Input, Dense, Dropout, Concatenate, Flatten, Conv2D, MaxPool2D
from keras.models import Model, Sequential
from keras.utils.np_utils import to_categorical
import matplotlib.pyplot as plt
import numpy as np
import glob
from PIL import Image
from keras.utils import np_utils
from sklearn.model_selection import train_test_split
from keras.preprocessing.image import load_img, img_to_array
classes = [
           "Dog", "Cat", "Rabbit"
           ]
X = []
Y = []
for index, classlabel in enumerate (classes):
    dir = "./images/" + classlabel
    files = glob.glob (dir + "/*.jpeg")
    for i, file in enumerate (files):
        image = Image.open (file)
        image = image.convert ("RGB")
        image = image.resize ((64, 64))
        data = np.asarray (image)
        X.append (data)
        Y.append (index)
X = np.array (X)
Y = np.array (Y)
X = X.astype ('float32') /255.0
# ★★ Fixed pre-processing from here
#Number of categories
num_classes = Y.max () --Y.min () + 1
# Divide image X into x1 and x2 and define a new correct label y
split_size = len (X) // 2
x1, x2, y1, y2 = train_test_split (X, Y, test_size = split_size, train_size = split_size)
# y is made from y1 and y2, which is the average value here. Is the average of dogs and rabbits cats?
y = (y1 + y2) // 2
x_train1, x_test1, x_train2, x_test2, y_train, y_test = train_test_split (x1, x2, y, test_size = 0.30)
# ★★ Corrected so far
Convert to # one-hot representation.
y_train = to_categorical (y_train, num_classes)
y_test = to_categorical (y_test, num_classes)
# ★★ Modify the model from here for the image
#Create a model.
model1 = Sequential ([[
    Conv2D (64,3, input_shape = (64,64,3)),
    Activation ('relu'),
    Conv2D (64,3),
    Activation ('relu'),
    MaxPool2D (pool_size = (2,2)),
    Flatten (),
    Dense (1024),
    Activation ('relu')
])
model2 = Sequential ([[
    Conv2D (64,3, input_shape = (64,64,3)),
    Activation ('relu'),
    Conv2D (64,3),
    Activation ('relu'),
    MaxPool2D (pool_size = (2,2)),
    Flatten (),
    Dense (1024),
    Activation ('relu')
])
concatenated = Concatenate (axis = -1) ([model1.output, model2.output])
# Tightly coupled
out = Dense (num_classes, activation ='softmax') (concatenated)
model = Model ([model1.input, model2.input], out)
# ★★ Corrected so far
#Model compilation
model.compile (loss ='binary_crossentropy', optimizer ='adam', metrics = ['acc'])
model.summary ()
# learn.
history = model.fit ([x_train1, x_train2], y_train, epochs = 40)