Several issues, the firstest is that you need the element-wise multiplication operator for P2, not matrix multiplication. These are the "dot" operators in Matlab...
P2=pi*d*r*(1+B*(S-S0)/d).*(Le-S+S0);
After that, with
Le=0.46;
S0=14.6;
then
S = S0:0.1:Le; --> S= 14.6: 0.1: 0.46;
which results in 0-length S since the end value is greater than the start.
There's something funky going on here w/ S as shown in your figure and how the equations are written it seems. It looks like the plot S is a normalized value 0:1 whereas you're trying to use some other length(?). Not sure how to grapple with that. If I simply use the figure normalized breakpoint that looks to correlate with your S0/100, I can get a plot of the shapes but the normalization is off...you'll have to figure out how that was actually done, precisely...
>> E=60000; d=0.014; n=4500; G=6.0; Le=0.46; B=0.05; r=3.0; S0=14.6;
>> S0=S0/100; % put back to scale of figure, roughly...
>> S=linspace(0,S0,20); % some points between origin and S0
>> P1 = sqrt(((pi^2*r*E*d^3*(1+n)*S)/2)+((pi^2*G*E*d^3)/2));
>> S2=linspace(S0,1,50); % points between S0 and 1
>> P2=pi*d*r*(1+B*(S2-S0)/d).*(Le-S2+S0);
>> hold all
>> plot(S2,P2*100+20) % arbitrary adjust P2 to get roughly same scaling...
results in
As noted, you'll need to figure out the actual scaling between the presented figure and the equations...there's some sleight-of-hand going on there it seems.
