syms f d m1 m2 G
m1 = 2e13;
m2 = 3e13;
G = 6.67e-11;
d1 = 0.01:0.001:0.050;
len = length(d1);
F = ones(len, 1);
tic
S = solve(f == G*m1*m2/(d*d), f);
for k = 1:len
F(k) = vpa(subs(S, d, d1(k)));
end
toc
d = 0.01:0.001:0.050;
figure
plot(d, F)
xlabel('d')
ylabel('F')
title('F~d')
legend ('d = 0.01:0.001:0.05')
tic
F1 = ones(len, 1);
for d = 0.01:0.001:0.05
F1(k) = 6.67e-11*2e13*3e13/(d*d);
end
toc