Plot function solution in 3D

I really would like to get some help draw a solution to this problem. I know c, x_0, sigma and sigma_H, i.e. c', delta and x_T is the variables to find.

Furthermore 0<c'<4, 1<delta<3 and 0.8<x_T<1

In addition I need to have x_T < f(c,delta)

I have tried

[c,delta,xT]= meshgrid(linspace(0,4,100),linspace(1,3,100),linspace(0.8,1,100));
V = straight.S0Func(optimalStraightCoupon,1, 1, straight.sigmaH)- straight.S0Func(c, delta, xT, straight.sigmaL); 
p = patch(isosurface(delta,xT,c,V,-0));
isonormals(c,delta,xT,V,p);
set(p,'FaceColor','red','EdgeColor','none');  
view(3);

Hope someone can help as I have been working with this problem for a couple of months and simply can not get the graph I am supposed to

This is what it should look like

further info

optimalStraightCoupon = fminbnd(@(c) -(straight.S0Func(c,1,1,straight.sigmaH)+straight.D0Func(c,1,1,straight.sigmaH)), 0, 4);

The straight class:

classdef Straight        
    properties
        sigmaL=0.2;
        sigmaH=0.3;
        x0=1;
        r=0.07;
        mu=0.05;
        alpha=0.15;
        tau=0.35;                  
        betaL;  
        betaH;
    end

      methods
          function beta = betaFunc(obj, sigma)
              beta = -(obj.mu-(1/2)*sigma^2+sqrt(2*obj.r*sigma^2+(obj.mu-(1/2)*sigma^2)^2))/sigma^2;
          end

          %Set both betas
          function obj = setBeta(obj)
              obj.betaL = obj.betaFunc(obj.sigmaL);
              obj.betaH = obj.betaFunc(obj.sigmaH);
          end

          %Calculate the default value of the firm at the cash level x         
          function L = LFunc(obj, x)
              L = (1-obj.alpha)*(1-obj.tau)*x/(obj.r-obj.mu);
          end   

          %Optimal default boundary given coupon, step-up and risk parameters          
          function xb = xbFunc(obj, c, delta, sigma) 
              beta = obj.betaFunc(sigma);
              xb = delta.*c.*(obj.r-obj.mu).*beta./(obj.r.*(beta-1));
          end

          %Optimal coupon for a straight bond with the high risk.
          function cPrime = cPrime(obj)
              xbHat = obj.xbFunc(1,1,obj.sigmaH);
              L = obj.LFunc(xbHat);
              cPrime = (((obj.tau)./(obj.r.*(1-obj.betaH))).*(((1-obj.tau).*((xbHat./(obj.r-obj.mu))-(1./obj.r))-(L-(1./obj.r))).^(-1)).*((obj.x0./xbHat).^(-obj.betaH))).^(-1./obj.betaH);
          end

          %Calculate S0 value     
          function S0 = S0Func(obj, c, delta, xT, sigma)
              beta = obj.betaFunc(sigma);
              xb = obj.xbFunc(c, delta, sigma);            
              S0 = (1-obj.tau).*(obj.x0./(obj.r-obj.mu)-c./obj.r-(xb./(obj.r-obj.mu)-c./obj.r).*((obj.x0./xb).^beta)-((delta-1).*c./obj.r).*((obj.x0./xT).^beta-((obj.x0./xb).^beta)));
          end         

          %Calculate D0 value
          function D0 = D0Func(obj, c, delta, xT, sigma)            
              beta = obj.betaFunc(sigma);
              xb = obj.xbFunc(c, delta, sigma);            
              L = obj.LFunc(xb);            
              D0 = c./obj.r+(L-c./obj.r).*((obj.x0./xb).^beta)+((delta-1).*c./obj.r).*((obj.x0./xT).^beta-((obj.x0./xb).^beta));
          end

Hope this give some sense of meaning, else do not hesitate to write.