sol = solve('h+r*x-x^3','x');

lw = 3;
%return

while 1
	h0 = input('Input value of h:');
	if isempty(h0)
        return
	end
    hold off
    ezsurf(sol(1),100)
	hidden off
	hold on
	
	ezsurf(sol(2),100)
	ezsurf(sol(3),100)
	axis([-6 6 -6 6 -4 4])
        alpha(0.5)

	if h0 == 0
        T = [0:0.001:6];
        plot3(zeros(size(T)),-T,zeros(size(T)),'LineWidth',lw); 
        plot3(zeros(size(T)),T,zeros(size(T)),'LineWidth',lw,'LineStyle','--');
        T = [-sqrt(6):0.02:sqrt(6)];
        plot3(zeros(size(T)),T.^2,T,'LineWidth',lw);
	else

        T = [-6:0.01:6];
        X = subs(sol(1),{'h','r'},{h0*ones(size(T)),T});
        ind = find(~imag(X));
        plot3(h0*ones(size(T(ind))),T(ind),X(ind),'LineWidth',lw);
        T = [-6:0.01:6];
        X = subs(sol(2),{'h','r'},{h0*ones(size(T)),T});
        ind = find(~imag(X));
        plot3(h0*ones(size(T(ind))),T(ind),X(ind),'LineWidth',lw);
        T = [-6:0.01:6];
        X = subs(sol(3),{'h','r'},{h0*ones(size(T)),T});
        ind = find(~imag(X));
        plot3(h0*ones(size(T(ind))),T(ind),X(ind),'LineWidth',lw,'LineStyle','--');
	end
    title(sprintf('h = %f',h0));
    view(-130,30);
end