for alp=2
    for B=0:0.1: 6
        f=@(theta)(B*(12-12*theta)*((4*alp-1)^2)-theta*(24*alp*(1-theta)+25*alp)*(12*(4*alp-1)-(2*alp-1)*(24-12*theta)+alp));
        m1=fzero(f,[0 1]);%(DU)

       g=@(theta)(B*12*(1-theta)*((4*alp-1)^2)-theta*(24+2*alp*(1-theta))*(24*alp-alp*(2*alp-1)*(1-theta)));
       m2=fzero(g,[0,1]);%(UD)

       h=@(theta)(B*12*(1-theta)*((4*alp-1)^2)-theta*(24*alp+25*alp*(1-theta))*(12*(4*alp-1)-(2*alp-1)*(24-12*theta)+alp*(1-theta))-...
           theta*(12+12*(1-theta)+2*alp*(1-theta))*(12*alp*(1-theta)-alp*(2*alp-1)*(1-theta)+12*alp));
       m3=fzero(h,[0 1]);%(DD)

       Y=max(12/(2*alp+12)*m1,m3);
       Z=max(2*alp/(2*alp+12)*m2,Y);
       if Z==12/(2*alp+12)*m1 && Z==2*alp/(2*alp+12)*m2
           pUU=1-(24*(alp-1)+2*alp*(alp-1))/(12*(4*alp-1));
            subplot(1,2,2),plot(B,pUU,'ko')%(UU)
            hold on
       else if Z==12/(2*alp+12)*m1 %(DU)
                pDU=1-(12*(1-m1)*(alp-1)+2*alp*(alp-1)+12*(alp-1))/(12*(4*alp-1));
            subplot(1,2,2), plot(B,pDU,'rp','LineWidth',0.1)
             hold on
           else if Z==2*alp/(2*alp+12)*m2 % (UD)
                   pUD=1-(24*(alp-1)+2*alp*(alp-1)*(1-m2))/(12*(4*alp-1));
              subplot(1,2,2),   plot(B,pUD,'g^','LineWidth',0.1)
                 hold on
               else pDD=1-((1-m3)*(alp-1)*(2*alp+12)+12*(alp-1))/(12*(4*alp-1)); % (DD)
               subplot(1,2,2),   plot(B,pDD,'*','LineWidth',0.1)
                  hold on
               end
           end  
       end
    end
end
legend('(U,U)','(D,U)','(D,D)')
xlabel('Budget B')
ylabel('Number of Patients (\Lambda)')