function [robots] = odesolve(robots,leader, t, numrobots)
m = 1;
C_0 = 1;
C_c = 0.7;
C_a = 50;
C_r = 90;
l_r = 12;
l_a = 2;
alpha = 0;
beta = 0.5;
maxtime = 200;
delta_t = 0.0001;
v(1:numrobots, 1, 1:2) = [robots(1:numrobots,t-1,3), robots(1:numrobots,t-1,4)];
for k = 1:maxtime
    for i = 1:numrobots
        term1 = [0,0];
        term2 = [0,0];
        term3 = [0,0];
        term4 = [0,0];
        U1 = 0;
        U2 = 0;
        for j = 1:(numrobots + 1)
            if j~=i
                if (j == numrobots + 1)
                    posx = leader(t-1,1);
                    posy = leader(t-1,2);
                    vx = leader(t-1,3);
                    vy = leader(t-1,4);
                else
                    posx = robots(j,t-1,1);
                    posy = robots(j,t-1,2);
                    vx = v(j,k,1);
                    vy = v(j,k,2);
                end
                n = -norm([robots(i,t-1,1) - posx, robots(i,t-1,2) - posy]);
                if (robots(i,t-1,1) > posx)
                    U1 = U1 + (-C_r/l_r)*exp(n/l_r) + (C_a/l_a)*exp(n/l_a);
                else
                    U1 = U1 + (C_r/l_r)*exp(n/l_r) - (C_a/l_a)*exp(n/l_a);
                end
                if (robots(i,t-1,2) > posy)
                    U2 = U2 + (-C_r/l_r)*exp(n/l_r) + (C_a/l_a)*exp(n/l_a);
                else
                    U2 = U2 + (C_r/l_r)*exp(n/l_r) - (C_a/l_a)*exp(n/l_a);
                end
                term2(1) = term2(1) + (1/m)*C_0*(vx - v(i,k,1));
                term2(2) = term2(2) + (1/m)*C_0*(vy - v(i,k,2));
            end
        end
        term1 = [-(1/m)*U1, -(1/m)*U2];
        term3(1) = (1/m)*(alpha - beta*(v(i,k,1)^2 + v(i,k,2)^2))*v(i,k,1);
        term3(2) = (1/m)*(alpha - beta*(v(i,k,1)^2 + v(i,k,2)^2))*v(i,k,2);
        leadvec = [leader(t-1,1) - robots(i,t-1,1), leader(t-1,2) - robots(i,t-1,2)];
        term4(1) = (1/m)*C_c*leadvec(1);
        term4(2) = (1/m)*C_c*leadvec(2);
        v(i,k+1,1) = v(i,k,1) + delta_t*(term1(1) + term2(1) + term3(1) + term4(1));
        v(i,k+1,2) = v(i,k,2) + delta_t*(term1(2) + term2(2) + term3(2) + term4(2));
    end
end

for i=1:numrobots
    robots(i,t,3) = v(i,end,1);
    robots(i,t,4) = v(i,end,2);
end

