@Preben Alsholm 

Thanks for your time, I will try some other tool then.

If anyone else has a walk around how to plot what I want in maple 2015 feel free to post it here.

@Preben Alsholm 

EDIT: I have license only for maple 2015. Any walkarounds for this issue?

I copy pasted your equations resulting in folowing worksheet:
restart; with(plots); f__r := proc (sigma__x, sigma__y, tau__xy, f__tx, f__ty, alpha) options operator, arrow; (1/2)*sigma__x-(1/2)*f__tx+(1/2)*sigma__y-(1/2)*f__ty+sqrt(((1/2)*sigma__x-(1/2)*f__tx-(1/2)*sigma__y+(1/2)*f__ty)^2+alpha*tau__xy^2) end proc; f__h := proc (sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u) options operator, arrow; sqrt((sigma__x*sigma__cy/sigma__cx+(1/2)*sigma__y*(1/f__45^2-1/f__cx^2-1/f__cy^2)*f__cx*f__cy)*sigma__x+((1/2)*sigma__x*(1/f__45^2-1/f__cx^2-1/f__cy^2)*f__cx*f__cy+sigma__y*sigma__cx/sigma__cy)*sigma__y+tau__xy^2*f__cx*f__cy/tau__u^2)-sqrt(sigma__cx*sigma__cy) end proc; yield := proc (sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u, f__tx, f__ty, alpha) if evalf(f__r(sigma__x, sigma__y, tau__xy, f__tx, f__ty, alpha)) < 0 then return f__h(sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u) else undefined end if end proc; plots:-implicitplot3d('yield(sigma__x, sigma__y, tau__xy, `$`(5, 9))' = 0, sigma__x = -10 .. 10, sigma__y = -10 .. 10, tau__xy = 0 .. 10, style = surfacecontour, numpoints = 100000, axes = normal);

and got:


any ideas?

@Preben Alsholm 
Sure this is whole worksheet:

restart: with(plots):

f__r := proc (sigma__x, sigma__y, tau__xy, f__tx, f__ty, alpha) options operator, arrow; (1/2)*sigma__x-(1/2)*f__tx+(1/2)*sigma__y-(1/2)*f__ty+sqrt(((1/2)*sigma__x-(1/2)*f__tx-(1/2)*sigma__y+(1/2)*f__ty)^2+alpha*tau__xy^2) end proc;

f__h := proc (sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u) options operator, arrow; sqrt((sigma__x*sigma__cy/sigma__cx+(1/2)*sigma__y*(1/f__45^2-1/f__cx^2-1/f__cy^2)*f__cx*f__cy)*sigma__x+((1/2)*sigma__x*(1/f__45^2-1/f__cx^2-1/f__cy^2)*f__cx*f__cy+sigma__y*sigma__cx/sigma__cy)*sigma__y+tau__xy^2*f__cx*f__cy/tau__u^2)-sqrt(sigma__cx*sigma__cy) end proc;

yield := proc (sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u, f__tx, f__ty, alpha) if eval(evalf(f__r(sigma__x, sigma__y, tau__xy, f__tx, f__ty, alpha)) < 0) then return f__h(sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u) else 9999999 end if end proc;

Explore(implicitplot3d('yield(sigma__x, sigma__y, tau__xy, sigma__cx, sigma__cy, f__45, f__cx, f__cy, tau__u, f__tx, f__ty, alpha)' = 0, sigma__x = -10 .. 10, sigma__y = -10 .. 10, tau__xy = 0 .. 10, style = surfacecontour, numpoints = 100000, axes = normal), parameters = [f__cx = 0.1e-3 .. 10, f__cy = 0.1e-3 .. 10, f__45 = 0.1e-3 .. 10, sigma__cx = 0.1e-3 .. 10, sigma__cy = 0.1e-3 .. 10, tau__u = 0.1e-3 .. 10, f__tx = 0.1e-3 .. 10, f__ty = 0.1e-3 .. 10, alpha = 0.1e-3 .. 10]);

EDIT:
Regarding two surfaces:

I have a cone and elipsoid:

These two surfaces should be displayed acording folowing conditions:

Display ellipsoid if function of cone is <0
Display cone if function of elipsoid is <0

In the end it should look similar to this (this is what my worksheet already produces):

But I need the folowing:

The new surface should be smooth [very important]
I would like to have different colors for different parts (elipsoid,cone) [if this is not possible it's not that important]

@Preben Alsholm 

Thanks, but using undifined I get an empty plot...

Here are my f__r and f__h, hope you can help.

Edit: Ideally I would like to have both surfaces ploted together exactly like in the image I showed in the initial post. though they both shoud be smooth as I later plot points and to see how these surfaces fit them.