clear all
h = 10;
h1 = h/2;
hzero= 0.3*h;

ctop=7;
c2=0;
for h2=0.9*hzero:0.02*hzero:1.1*hzero
     for c1=0.9*ctop:0.02*ctop:1.1*ctop
           a=(c1-c2)/(h1^2-h2^2);
          Tmax=c1-a*h1^2;
           y_co = h/2:(-h/40):-h/2;
          sigma = a.*y_co.^2 + Tmax;
          plot(sigma, y_co);
          hold on
     end 
end