Este manuscrito presenta una técnica discreta para estimar la solución de un condensado de Bose-Einstein de dos componentes doble fraccional. El sistema consta de dos ecuaciones diferenciales parciales parabólicas no lineales acopladas cuyas soluciones son dos funciones complejas, y las derivadas fraccionarias espaciales se interpretan en el sentido de Riesz. Se imponen datos iniciales y de contorno de Dirichlet homogéneos en un dominio espacial multidimensional. Para aproximar las soluciones, empleamos una metodología de diferencias finitas. Establecemos rigurosamente la existencia de soluciones numéricas junto con las principales propiedades numéricas. Concretamente, demostramos que el esquema es consistente tanto en espacio como en tiempo, además de ser estable y convergente. Se presentan simulaciones numéricas en el escenario unidimensional para mostrar el rendimiento del esquema. Por conveniencia, se proporciona un código MATLAB del modelo numérico en el apéndice al final de este trabajo.