function poly()
% Polynom 3. Grades
wdh = 'j';
while wdh == 'j'
	clc;
	disp('****  Polynom dritten Grades: y=ax³+bx²+cx+d  ****');
	disp(' ');
	disp('(1)	Beispiel 1:  a=-1 b=3 c=0 d=-1');
	disp('(2)	Beispiel 2:  a=3 b=-4 c=-4 d=0');
	disp('(3)	Koeffizienten eingeben:  a b c d');
	disp('(4)	Programm beenden');
	disp('');
	switch input('Auswahl: ')
	case 1, clc
        a=-1;, b=3;, c=0;, d=-1;
        poly3_rechnung(a,b,c,d)      
        disp(' ');
        disp('Weiter mit Enter');
        pause;
        close;
	case 2, clc
        a=3;, b=-4;, c=-4;, d=0;
        poly3_rechnung(a,b,c,d)       
        disp(' ');
        disp('Weiter mit Enter');
        pause;
        close;
	case 3, clc, disp('Koeffizienten eingeben');
        a=input('a: ');, b=input('b: ');, c=input('c: ');, d=input('d: ');
        clc;
        poly3_rechnung(a,b,c,d)
        disp(' ');
        disp('Weiter mit Enter');
        pause;
        close;
	case 4, wdh = 'n';
	otherwise disp('');
	     disp('Wählen Sie eine gültige Nummer !');
	     disp('Weiter mit "Enter"!');
         pause;
         close;
	end;	
end;

function poly3_rechnung(a,b,c,d)
% wenn a=0 dann nich gut
if a==0
    disp('a=0 => kein Polynom dritten Grades => Berechnung nicht möglich');
else
    % Extrempunkte: 1,2 : 1. Ableitung der Funktion = 0 => quadratische Gleichung lösen 
    exp2=-(2/3*b/a/2)-sqrt((2/3*b/a/2)^2-1/3*c/a);      
    exp1=-(2/3*b/a/2)+sqrt((2/3*b/a/2)^2-1/3*c/a);      
    % Wendepunkt: 2. Ableitung = 0 => nach x umgestellt 
    W=-(2/6*b/a);
    % Festlegung Extrempunkte Hoch- bzw Tiefpunkt
    if (6*a*exp1+2*b)<0                              
        H=exp1;
        T=exp2;
    else
        T=exp1;
        H=exp2;
    end
    % y-Werte für Extrempunkte und Wendepunkt
    yH=a*H^3+b*H^2+c*H+d;
    yT=a*T^3+b*T^2+c*T+d;
    yW=a*W^3+b*W^2+c*W+d;
    % Vorzeichenvariable für a
    if a>0
        vz=1;
    else
        vz=-1;
    end
    % Wenn ein Extrempunkt imaginäre Zahl dann nur eine Nullstelle ( unter der Wurzel kleiner 0 )
    if ((2/3*b/a/2)^2-1/3*c/a)<=0 
        if yW==0
            nst1=0;
            nst2=[];
            nst3=[];
            H=[];
            T=[];
            yH=[];
            yT=[];
        else        
            nst1=poly3_iteration(a,b,c,d,W+1);
            nst2=[];
            nst3=[];
            H=[];
            T=[];
            yH=[];
            yT=[];
        end
        i=abs(W-nst1);
        j=abs(yW);
        if W==nst1
            i=abs(10*a);
            j=10*abs(a*(nst1+i/10).^3+b*(nst1+i/10).^2+c*(nst1+i/10)+d);
        end
        if a>=0
            if yW>0
                x=nst1-i/10:i/100:W+i/10;
                ymin=0-j/10;
                ymax=yW+j/10;
            else
                x=W-i/10:i/100:nst1+i/10;
                ymin=yW-j/10;
                ymax=0+j/10;
            end   
        else
            if yW>0
                x=W-i/10:i/100:nst1+i/10;
                ymin=0-j/10;
                ymax=yW+j/10;
            else
                x=nst1-i/10:i/100:W+i/10;
                ymin=yW-j/10;
                ymax=0+j/10;
            end              
        end
    % beide Extrempunkte über X-Achse
    elseif yT>0 
        nst1=poly3_iteration(a,b,c,d,H-vz*abs(T-H)/2);
        nst2=[];
        nst3=[];
        i=abs(T-nst1);
        j=abs(yH);
        ymin=0-j/10;
        ymax=yH+j/10;
        if a>=0
            x=nst1-i/10:i/100:T+i/10;
        else
            x=T-i/10:i/100:nst1+i/10;
        end
    % beide Extrempunkte unter X-Achse
    elseif yH<0 
        nst1=poly3_iteration(a,b,c,d,T+vz*abs(T-H)/2);
        nst2=[];
        nst3=[];
        i=abs(H-nst1);
        j=abs(yT);
        ymin=yT-j/10;
        ymax=0+j/10;
        if a>=0
            x=H-i/10:i/100:nst1+i/10;      
        else
            x=nst1-i/10:i/100:H+i/10;
        end
    % Extrempunkte und Wendepunkt in einem Punkt (0,0)
    elseif H==T 
        if yH==0
            nst1=H;
            nst2=[];
            nst3=[];
            x=nst1-1:0.01:nst1+1;
            ymin=0-1;
            ymax=0+1;
        end
    % Hochpunkt = Nullstelle    
    elseif yH==0
        nst1=poly3_iteration(a,b,c,d,T+vz*abs(T-H)/2);
        nst2=H;
        nst3=[];
        i=abs(H-nst1);
        j=abs(yT);
        ymin=yT-j/10;
        ymax=0+j/10;
        if a>=0
            x=H-i/10:i/100:nst1+i/10;
        else
            x=nst1-i/10:i/100:H+i/10;
        end
    % Tiefpunkt = Nullstelle    
    elseif yT==0
        nst1=poly3_iteration(a,b,c,d,H-vz*abs(T-H)/2);
        nst2=T;
        nst3=[];
        i=abs(T-nst1);
        j=abs(yH);
        ymin=0-j/10;
        ymax=yH+j/10;
        if a>=0
            x=nst1-i/10:i/100:T+i/10;
        else
            x=T-i/10:i/100:nst1+i/10;
        end
    % 3 Nullstellen  
    else 
        nst1=poly3_iteration(a,b,c,d,H-vz*abs(T-H)/2);     
        if yW==0
            nst2=0;
        else
            nst2=poly3_iteration(a,b,c,d,W);
        end
        nst3=poly3_iteration(a,b,c,d,T+vz*abs(T-H)/2);
        i=abs(nst3-nst1);
        j=abs(yH-yT);
        ymin=yT-j/10;
        ymax=yH+j/10;
        if a>=0
            x=nst1-i/10:i/100:nst3+i/10;
        else
            x=nst3-i/10:i/100:nst1+i/10;
        end
    end 
    % y-Werte der Nullstellen nur um zu prüfen
    ynst1=a*nst1^3+b*nst1^2+c*nst1+d;
    ynst2=a*nst2^3+b*nst2^2+c*nst2+d;
    ynst3=a*nst3^3+b*nst3^2+c*nst3+d;
    % alle y_Werte für zugehörige x-Werte
    y=a*x.^3+b*x.^2+c*x+d;
    % Schreiben der Funktion
    if a==1
        va='x³';
    elseif a==-1
        va='-x³';
    elseif a>0
        va=[num2str(a),'x³'];
    elseif b<0
        va=[num2str(a),'x³'];
    else
        va='';
    end
    if b==1
        vb='+x²';
    elseif b==-1
        vb='-x²';
    elseif b>0
        vb=['+',num2str(b),'x²'];
    elseif b<0
        vb=[num2str(b),'x²'];
    else
        vb='';
    end
    if c==1
        vc='+x';
    elseif c==-1
        vc='-x';
    elseif c>0
        vc=['+',num2str(c),'x'];
    elseif c<0
        vc=[num2str(c),'x'];
    else
        vc='';
    end
    if d>0
        vd=['+',num2str(d)];
    elseif d<0
        vd=num2str(d);
    else
        vd='';
    end
    disp(' ');
    disp(['y=f(x)=',va,vb,vc,vd]);
    disp(' ');
    disp(['X-Null 1   :   ',num2str(nst1)]);
    disp(['X-Null 2   :   ',num2str(nst2)]);
    disp(['X-Null 3   :   ',num2str(nst3)]);
    % Sattelpunkt wenn W in erste Ableitung eingesetzt 0 ergiebt    
    if (3*a*W^2+2*b*W+c)==0
        disp(['Sattelpunkt  :   (',num2str(W),' , ',num2str(yW),')']);
    else
        disp(['Hochpunkt  :   (',num2str(H),' , ',num2str(yH),')']);
        disp(['Tiefpunkt  :   (',num2str(T),' , ',num2str(yT),')']);
        disp(['Wendepunkt :   (',num2str(W),' , ',num2str(yW),')']);
    end
    % grafische Ausgabe der Funktion
    plot(x,x*0,'k',x,y,nst1,ynst1,'ro',W,yW,'go',H,yH,'mo',T,yT,'ko',nst2,ynst2,'ro',nst3,ynst3,'ro'); 
    set(1,'Position',[0 50 600 400])
    grid on
    axis([x(1) x(end) ymin ymax])
    legend('x-Achse','Graph','Nullstellen','Wendepunkt','Hochpunkt','Tiefpunkt',0)
    title(['y=f(x)=',va,vb,vc,vd])
    xlabel('x-Achse')
    ylabel('y-Achse')
end

function y=poly3_iteration(a,b,c,d,x0)
% Funktion zur iterativen Nullstellenberechnung
% eines Polynoms 3.ten Grades ax^3+bx^2+cx+d = 0 
% nach dem Newton - Verfahren bei gegebenen Startwert x0.
% Aufruf: y=poly3_iteration(a,b,c,d,x0).
% Bei erfolgreicher Iteration liefert y eine Nullstelle,
% andernfalls ist y eine leere Matrix.
	anz = 0;   
	f = a*x0^3+b*x0^2+c*x0+d;
    while abs(f) > 1e-8	% Iterationsschleife : Anfang
		f1 = 3*a*x0^2+2*b*x0 + c;
		x0 = x0 - f/f1;
		anz = anz + 1;
		f = a*x0^3+b*x0^2+c*x0+d;
        if anz>30
            break
        end
    end
		y=x0;
end % Iterationsschleife : Ende