function bead

clf reset;
hold on;

x=-5:.01:5;

for E=5:5:50
    
    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

hold on;
subplot(2,1,1);
xlabel('x');ylabel('p');
hold on;
subplot(2,1,2);
xlabel('x');ylabel('x_dot');