% original solution
alpha=3.5;
[x,y]=rk3(@pendulum_ode,[0,10],[0;alpha],500);
figure(1)
plot(x,y(1,:),'k','linewidth',2)
print -dpng pendulum2a.png
close(1)

% solution with small alpha
alpha=0.035;
[x,y]=rk3(@pendulum_ode,[0,10],[0;alpha],500);
% exact solution
exact=alpha/sqrt(3)*sin(sqrt(3)*x);
figure(2)
plot(x,exact,'g','linewidth',8);
hold on
plot(x,y(1,:),'k','linewidth',2);
hold off
title('alpha = 0.035')
print -dpng pendulum2b.png

% solution with medium alpha
alpha=0.35;
[x,y]=rk3(@pendulum_ode,[0,10],[0;alpha],500);
% exact solution
exact=alpha/sqrt(3)*sin(sqrt(3)*x);
figure(3)
plot(x,exact,'g','linewidth',8);
hold on
plot(x,y(1,:),'k','linewidth',2);
hold off
title('alpha = 0.35')
print -dpng pendulum2c.png

% solution with larger alpha
alpha=2;
[x,y]=rk3(@pendulum_ode,[0,10],[0;alpha],500);
% exact solution
exact=alpha/sqrt(3)*sin(sqrt(3)*x);
figure(4)
plot(x,exact,'g','linewidth',8);
hold on
plot(x,y(1,:),'k','linewidth',2);
hold off
title('alpha = 2')
print -dpng pendulum2d.png

% solution with larger yet alpha
alpha=3.2;
[x,y]=rk3(@pendulum_ode,[0,10],[0;alpha],500);
% exact solution
exact=alpha/sqrt(3)*sin(sqrt(3)*x);
figure(5)
plot(x,exact,'g','linewidth',8);
hold on
plot(x,y(1,:),'k','linewidth',2);
hold off
title('alpha=3.2')
print -dpng pendulum2e.png

% another nonlinear solution
alpha=3.6;
[x,y]=rk3(@pendulum_ode,[0,10],[0;alpha],500);
[xx,yy]=rk3(@(x,y) [y(2);-3*(y(1)-y(1)^3/6+y(1)^5/120-y(1)^7/5040)],[0,2.5],[0;alpha],1000);
figure(6)
plot(xx,yy(1,:),'k','linewidth',2);
hold on
plot(x(1:150),y(1,(1:150)),'r','linewidth',2);
hold off
legend('similar equation','original')
title('similar nonlinear equation')
print -dpng pendulum2f.png
