function myarea=ngon(sides,shapes,n)

R=1;G=0;B=0;
color=[R G B];

Q=[1];
P=[0];
theta=2*pi/sides;

for j=2:sides
	Q=[Q;cos( (j - 1) * theta )];
	P=[P;sin( (j - 1) * theta )];
end

myarea=polyarea(Q,P);

figure;
hold on;
fill(Q,P,color);
axis equal;
axis off;

h=0.01;
%n=2
N=5;
if nargin < 2;shapes=20;end

qstring='Which solver would you like to use?';
ans1='Euler''s Method';
ans2='Euler Palais';
ans3='Runge-Kutta';
title='Integration Method?';

answer=questdlg(qstring,title,ans1,ans2,ans3,ans1);

for i=1:shapes;
	for j=1:sides
		q=Q(j);
		p=P(j);
		if strcmp(answer,ans1);
			for k=1:N
				qtemp=q;
				% Euler for Q
				q=q+h*p;
				% Euler for P
				p=p-h*qtemp^(2*n-1);
			end
		elseif strcmp(answer,ans2);
			for k=1:N
				% Euler for Q
				q=q+h*p;
				% Euler for P
				p=p-h*q^(2*n-1);
			end
		elseif strcmp(answer,ans3);
			for k=1:N

				qtemp=q;

				% rk4 for Q
				qh1 = h * p;
				qh2 = h * (p + (qh1) / 2);
				qh3 = h * (p + (qh2) / 2);
				qh4 = h * (p + qh3);
				q = q + (qh1) / 6 + (qh2) / 3 + (qh3) / 3 + (qh4) / 6;

				% rk4 for P
				ph1 = h * qtemp^(2*n-1);
				ph2 = h * (qtemp^(2*n-1) + (ph1) / 2);
				ph3 = h * (qtemp^(2*n-1) + (ph2) / 2);
				ph4 = h * (qtemp^(2*n-1) + ph3);
				p = p - (ph1) / 6 - (ph2) / 3 - (ph3) / 3 - (ph4) / 6;
			end
		end		
		Q(j)=q;
		P(j)=p;
	end
	myarea(i+1,:)=polyarea(Q,P);

%   SCALE FOR IDENTICALY OREINTED PLOT

%	scale=sqrt( myarea(i+1,:)/myarea(i) );
%	Qdraw=[scale * 1];
%	Pdraw=[0];
%	for j=2:sides
%		Qdraw=[Qdraw;scale * cos( (j - 1) * theta )];
%		Pdraw=[Pdraw;scale * sin( (j - 1) * theta )];
%	end
%	R=1-i/shapes;
%	fill(Qdraw,Pdraw,[1-i/shapes 0 0 ]);

	R=1-i/shapes;
	fill(Q,P,[R G B]);shg;

end

axis off;
hold off;
rotate3d on;
shg;