function bead

clf reset;
hold on;

x=-1:.01:1;

for E=.1:.5:10
    
    hold on;
    subplot(2,1,1);
    P_plus=sqrt(2*cosh(x).^2*E-9.8*(cosh(x)).^3);
    P_minus=-sqrt(2*cosh(x).^2*E-9.8*cosh(x).^3);
    plot(x,P_plus);
    plot(x,P_minus);
    
    hold on;
    subplot(2,1,2);
    x_dot_plus=sqrt(2*(E-9.8*cosh(x))/(cosh(x).^2));
    x_dot_minus=-sqrt(2*(E-9.8*cosh(x))/(cosh(x).^2));
    plot(x,x_dot_plus);
    plot(x,x_dot_minus);
    
end