#analytische und mumerische Lösung der DGL: dE/dt=-k1E+k2*E(N-E)
#k2=n*k1/N, n>1, Element von R , da k2 nur schwer bestimmbar ist
k1=5/73;
N=83e6;
a=1;
n=[1.5,2,3,5,10,100];
tmax=600;
delta=0.1;
X=0:delta:tmax;

function retval = wendeTangente(m,b,Emax,tt,tr,x)
  if(x >= tt && x <= tt+tr ) 
    retval=m*x+b;
  else 
     if(x <tt) 
      retval=0;
     else
       retval=Emax*x./x;
     endif;  
  endif;
endfunction;
  

erkrankte=@(N,a,k1,n,t) (n-1)*a*N./(n*a+(N*(n-1)-n*a)*exp(-k1*(n-1)*t));

Y=[];
for k = 1:length(n),
  y=@(X) erkrankte(N,a,k1,n(k),X);
  Y=[Y,y(X)'];
end; 



function retval = decorateSubplot(tmax)
   xlabel('t/Tage');
   ylabel('Erkrankte')
   xlim([0,tmax]);
   ylim([0,85e6]);
   set(gca,'XTick',0:25:tmax);
   set(gca,'YTick',0:10e6:90e6);
   grid on;
   retval=1;
endfunction;

subplot(3,1,1);
plot(X,Y);
decorateSubplot(tmax);
title("analytische Lösung");

Y=[];
twende=@(N,a,k1,n) -1/(k1*(n-1))*reallog(a*n/(N*(n-1) -n*a));
EPunkt=@(N,a,k1,n,t) N*a*k1*(n-1)^2*(N*(n-1)-n*a)*exp(-k1*(n-1)*t)/((n*a+(N*(n-1)-n*a)*exp(-k1*(n-1)*t))^2);
#fprintf("n;tw;E;limE;m;b;tt;tr\n");
for k = 1:length(n),
  tw=twende(N,a,k1,n(k));
  E=erkrankte(N,a,k1,n(k),tw);
  Emax=(n(k)-1)*N/n(k);
  m=EPunkt(N,a,k1,n(k),tw);
  b=E-m*tw;
  tt=-b/m;
  tr=Emax/m;
  #fprintf("n=%d, tw=%d, E=%d, limE=%d, m=%d, b=%d, tt=%d, tr=%d\n", n(k), tw, E, Emax, m, b, tt, tr );
  #fprintf("%d;%d;%d;%d;%d;%d;%d;%d\n", n(k), tw, E, Emax, m, b, tt, tr );
  for l=1:length(X),
     Y(k,l)=wendeTangente(m,b,Emax,tt,tr,X(l));
  end;
  # y=@(x) wendeTangente(m,b,Emax,tt,tr,x);
  # Y=[Y,y(X)'];
  #fprintf("%d %d\n" , length(y(X)'), length(X));
 
end; 

subplot(3,1,2);
plot(X,Y);
decorateSubplot(tmax);
title("Wendetangenten");

Y=[];
for k = 1:length(n),
  k2=k1*n(k)/N;
  erkranktePunkt=@(N,k1,k2,y,x) [-k1*y(1)+k2*y(1)*(N-y(1))];
  yPunkt=@(y,x) erkranktePunkt(N,k1,k2,y,x);
  Y=[Y, lsode(yPunkt,[a],X)];
end;  

subplot(3,1,3);
plot(X,Y);
decorateSubplot(tmax);
title("numerische Lösung lsode");