num_robots = 1000;
u = [1; 6];
X = zeros(3, num_robots);
X_prime = zeros(size(X));
for i = 1:num_robots
    X_prime(:, i) = motion_model(X(:, i), u);
end
X_prime

%

plot(X_prime(1,:), X_prime(2,:), 'b.')
axis equal

%

function x_prime = motion_model(x, u)
    dt = 0.1;
    a = [0.002 0.002 0.2 0.2];
    v = u(1);
    w = u(2);
    u_hat = u + [a(1) * v^2 + a(2) * w^2;...
        a(3) * v^2 + a(4) * w^2]...
    .* randn(2, 1);
    v = u_hat(1);
    w = u_hat(2);
    theta = x(3);
    x_prime = x +...
        [-v/w * sin(theta) + v/w * sin(theta + w * dt);...
        v/w * cos(theta) - v/w * cos(theta + w * dt);...
        w * dt];
end