clear;
E=zeros(1,100);
pi=3.1415926;
c=0.3;
u=zeros(1,101);
a=zeros(1,101);
d=zeros(1,100);
d2=zeros(1,100);
u(101)=0;a(101)=0;a(1)=1;
format long;
e=1403.765;
max=e;
min=e;

for x=100:-1:2
   u(x)=pi.*(2.*x-1)+(1-a(x+1)).*u(x+1);
   m=e-c.*(u(x)/(pi*(2*x-1))-1)-u(x)/(pi*(2*x-1));
   n=u(x)/(pi*(2*x-1))*(x.*x-1);
   a(x)=m./n;  
end
u(1)=pi.*(2.*x-1)+(1-a(2)).*u(2);
for x=1:100
   s(x)=a(x);
   d(x)=(1+u(x+1).*(1-a(x+1))./(pi.*(2.*x-1)));
   d2(x)=u(x+1).*(1-a(x+1))./(pi.*(2.*x-1));
   e1=c.*d2(x);
   e2=a(x).*x.*x.*d(x);
   e3=(1-a(x)).*d(x);
   E(x)=e1+e2+e3;
   if(max<E(x))
       max=E(x);
       i=x;
   end
   if(min>E(x))
       min=E(x);
       j=x;
   end
   
end

plot(s);
xlabel('D');
ylabel('a');
axis([0 120 0 1]);

%waitforbuttonpress;
figure;
format long;
E
max
min
z=max/min
i
j
plot(E);
xlabel('D');
ylabel('E');
axis([0 120 1403.76 1403.8]);