How do I run the following code in paralle

I want to make this code to run in parallel for each value of a, and b= 0:0.001:1, and there is a normalized condition that a^2+b^2+c^=1, so c=sqrt(1-a*a-b*b), in order to make plot 3D (fmax vs a, and b) . as this code run perfectly for just one value as you can see below. Any help would be appreciated. Thanks

m=0;
th1max=0;
th2max=0;
th3max=0;
th4max=0;
th5max=0;
th6max=0;
% 
step = 0.3;
for th1=(0:step:2)*pi
for th2=(0:step:2)*pi
for th3=(0:step:2)*pi
for th4=(0:step:2)*pi 
for th5=(0:step:2)*pi
for th6=(0:step:2)*pi    
a=.57;b=.57;c=sqrt(1-a*a-b*b);
p1=-sin(((th2-th1)/2)+((th4-th3)/2)+((th6-th5)/2))+sin(((th4-th3)/2)+((th6-th5)/2)-((th2-th1)/2))+sin(((th2-th1)/2)+((th6-th5)/2)-((th4-th3)/2))+sin(((th2-th1)/2)+((th4-th3)/2)-((th6-th5)/2));
p2=sin(((th2-th1)/2)+((th4-th3)/2)+((th6-th5)/2))+sin(((th4-th3)/2)+((th6-th5)/2)-((th2-th1)/2))-sin(((th2-th1)/2)+((th6-th5)/2)-((th4-th3)/2))+sin(((th2-th1)/2)+((th4-th3)/2)-((th6-th5)/2));
p3=sin(((th2-th1)/2)+((th4-th3)/2)+((th6-th5)/2))-sin(((th4-th3)/2)+((th6-th5)/2)-((th2-th1)/2))+sin(((th2-th1)/2)+((th6-th5)/2)-((th4-th3)/2))+sin(((th2-th1)/2)+((th4-th3)/2)-((th6-th5)/2));
p4=sin(((th2-th1)/2)+((th4-th3)/2)+((th6-th5)/2))+sin(((th4-th3)/2)+((th6-th5)/2)-((th2-th1)/2))+sin(((th2-th1)/2)+((th6-th5)/2)-((th4-th3)/2))-sin(((th2-th1)/2)+((th4-th3)/2)-((th6-th5)/2));
f=p1+2*a*c*p2+2*a*b*p3+2*b*c*p4;
if f>m
    m=f;
th1max=th1;
th2max=th2;
th3max=th3;
th4max=th4;
th5max=th5;
th6max=th6;
end
%display (f);
end 
end
end
end
end
end
m
th1max
th2max
th3max
th4max
th5max
th6max
**The result*** 

m =
      4.3530
th1max =
       0
th2max =
      1.8850
th3max =
       0
th4max =
      1.8850
th5max =
      1.8850
th6max =
      3.7699

[Walter Roberson]

Duplicates http://uk.mathworks.com/matlabcentral/answers/260783-how-to-run-the-following-code-in-parallel