how to prevent integration error in ode45 MATLAB

Question

i have a problem in my code for simulating the path of objects in a "4 body problem".

the goal is to solve differential equations for these 4 objects to obtain their position and velocity during time. the ub is a vector containing x,y,z,vx,vy,vz.

this is my MATLAB code:

clear;clc;close all;
ms1 = 10^30;ms2 = 10^30;mp1 = 10^6;mp2 = 10^6;
masses = [ms1;ms2;mp1;mp2];
G = 6.6742*10^11;
time = 30*86400;
%[r,v]
us1_0 = [10^10; 0; 0; 0; 40000; 0];
us2_0 = [-10^10; 0; 0; 0; -40000; 0];
up1_0 = [1.5*10^10; 0; 0; 0; 160000; 0];
up2_0 = [-1.5*10^10; 0; 0; 0; -160000; 0];
initials = [us1_0; us2_0; up1_0; up2_0];


[t, result] = ode45(@(t,ub) eq(t, ub, masses, G), [0,time], initials);

figure
x1 = result(:,1);y1 = result(:,2);z1 = result(:,3);
x2 = result(:,7);y2 = result(:,8);z2 = result(:,9);
x3 = result(:,13);y3 = result(:,14);z3 = result(:,15);
x4 = result(:,19);y4 = result(:,20);z4 = result(:,21);

grid on;
plot3(x1,y1,z1);hold on
plot3(x2,y2,z2);hold on
plot3(x3,y3,z3);hold on
plot3(x4,y4,z4);
legend('star1','star2','planet1','planet2')


function dubdt = eq(~, ub, masses, G)
    dubdt = zeros(24,1);
    dubdt(1:3, 1) = ub(4:6, 1);
    dubdt(7:9, 1) = ub(10:12, 1);
    dubdt(13:15, 1) = ub(16:18, 1);
    dubdt(19:21, 1) = ub(22:24, 1);
    for i = 1:4
        for j = 1:4
            if i ~= j 
                dubdt((6*i-2):(6*i),1) = dubdt((6*i-2):(6*i),1) + -G .* masses(j) .* (ub((6*i-5):(6*i-3), 1) - ub((6*j-5):(6*j-3), 1)) ...
                    /(norm(ub((6*i-5):(6*i-3), 1) - ub((6*j-5):(6*j-3), 1))).^3;
            end
         end
    end
end

when i try to run this, the error i mention below arises:

Warning: Failure at t=4.845555e-07.  Unable to meet integration tolerances without reducing the step size
below the smallest value allowed (1.694066e-21) at time t. 
> In ode45 (line 353)
In galaxy (line 14) 

what can i do to solve this?, thanks.