%DOPT  Optimal drag coefficient
%
C=zeros(5,9);
C(1,1)=1/12; C(2,1)=3/420; C(2,2)=-14/420;
C(2,3)=14/420; C(3,1)=1/1848; C(3,2)=-12/1848;
C(3,3)=45/1848; C(3,4)=-66/1848; C(3,5)=33/1848;
C(4,1)=1/25740; C(4,2)=-22/25740; C(4,3)=165/25740;
C(4,4)=-572/25740; C(4,5)=1001/25740; C(4,6)=-858/25740;
C(4,7)=286/25740; C(5,1)=3/1108536; C(5,2)=-104/1108536;
C(5,3)=1274/1108536; C(5,4)=-7644/1108536;
C(5,5)=25480/1108536; C(5,6)=-49504/1108536;
C(5,7)=55692/1108536; C(5,8)=-33592/1108536;
C(5,9)=8398/1108536;
for N=4:2:10
  P=P_matrix(N);
  A=P(3:N+1,3:N+1); b=2*P(1,3:N+1)';
  cond(A)
  c=A\b
  s1=sum(P(1,3:2:N+1)'.*c(1:2:N-1));
  s2=0;
  for k=3:2:N+1
    for l=3:2:N+1
      s2=s2+P(k,l)*c(k-2)*c(l-2);
    end
  end
  D=6-4*s1+s2
  hold on
  x=linspace(0,1)';
  S=x.*(1-x);
  for k=2:2:N
    p=zeros(size(x,1),1);
    for l=1:k-1
      p=p+C(k/2,l)*x.^(l-1);
    end
    S=S+c(k-1)*(x.*(1-x)).^2.*p;
  end 
  S=.1*S;
  R=sqrt(S/pi);
  if N==10
    plot(x,S,':')
    axis('equal')
    plot(x,R,'-')
    plot(x,-R,'-')
  end
end
