function [T,Y] = mzv(msintervall,startwerte,f,fy,r,dr_y0,dr_ytf,param)

% Anne Koehler, Andreas Sommer, Ilia Gherman
% 18.11.2004
%
% Mehrzielverfahren
% 
% 
%
% [T,Y] = mzv(msintervall,startwerte,f,fy,r,dr_y0,dr_ytf)
%
% msintervall   --> Gitter(-punkte)
% startwerte    --> Matrix, die Startwerte fuer s enthaelt  (als Spalten)
% f             --> Funktionsname fuer die "rechte Seite"
% fy            --> Funktionsname fuer die Ableitung der "rechten Seite" nach y
% r             --> Funktionsname fuer die Funktion der Randbedingungen
% dr_y0         --> Funktionsname fuer Ableitung der Randbedingungen nach y0
% dr_ytf        --> Funktionsname fuer Ableitung der Randbedingungen nach y(t_f)
% z = (y, G)    --> Loesungs-Vektor fuer y und G 
%
% T             --> Zeitpunkte der Loesung
% Y             --> Loesung der DGL

% Dimension der DGL und Anzahl der MS-Intervalle
[n,m]=size(startwerte);


% Multiple-Shooting
epsilon = 0.001;

delta_s=1+epsilon; % nur damit die while zum 1. mal ausgefuehrt wird

while (max(abs(delta_s)) > epsilon) 
   Y=[]; % hier speichern wir die gesamte Loesung ab
   T=[]; % 
   G=[]; % 
   % Single-Shooting auf jedem Teilintervall
   for i=1:m
       % Startwerte fuer zusammengesetzte DGL  
       z_null(1:n) = startwerte(:,i);
       z_null(n+1: n*n+n) = reshape(eye(n), n*n, 1);
       z_null=reshape(z_null,(n+1)*n,1);
            
       % Zusammengesetzte DGL loesen
       [T_temp,Z_temp] = ode15s(@h,[msintervall(i),msintervall(i+1)], z_null, [], f, fy, n,param);

       Y_temp = Z_temp(:,1:n);
       T = [T;T_temp];
       Y = [Y;Y_temp];
       G(i,:,:) = reshape(Z_temp(length(T_temp),(n+1):end)', n, n);
       y_ende(:,i) = Y_temp(length(T_temp),:)';
       
       % Ausgabe (nur sinnvoll, bei n=1 oder n=2, bzw. nur y1 plotten)
       plot(T_temp,Y_temp(:,1)) 
       hold on
       plot(T_temp, Y_temp(:,2), 'r');
       u = Y_temp(:,5)./(2.*Y_temp(:,5).*Y_temp(:,3));
       plot(T_temp, u, 'g');
       %disp('Taste...');
       %pause;
   end;
   pause(.3);
   clf;
   
   %Aufbau der MS-Matrix in Sparse-Form
   zeilen=[];
   spalten=[];
   werte=[];
   for i=1:m-1
        for j=1:n
            temp=n*(i-1);
            zeilen=[zeilen,[temp+1:1:temp+n]];
            spalten=[spalten,(temp+j)*ones(1,n)];
            werte=[werte,reshape(G(i,:,j),1,n)];
        end
        zeilen=[zeilen,[temp+1:1:temp+n]];
        spalten=[spalten,[n*i+1:1:n*i+n]];
        werte=[werte,-ones(1,n)];
    end

   for j=1:n
       zeilen=[zeilen,[n*(m-1)+1:1:n*m]];
       spalten=[spalten,j*ones(1,n)];
   end
   werte=[werte,reshape(feval(dr_y0,Y(1,:)',Y(end,:)',param),1,n*n)];
   for j=1:n
       zeilen=[zeilen,[n*(m-1)+1:1:n*m]];
       spalten=[spalten,(n*(m-1)+j)*ones(1,n)];
   end
   G_m = reshape(G(m,:,:),n,n);
   werte=[werte,reshape(feval(dr_ytf, Y(1,:)', Y(end,:)',param)*G_m,1,n*n)];
   A_ms_sparse=sparse(zeilen,spalten,werte,m*n,m*n);
              
   % Hier der 'alte' Aufbau der MS-Matrix,
   % nicht in Sparse-Variante
   %---------------------------------------------------
   % Aufbau der MS-Matrix (m Spalten und m Zeilen mit den (n x n) - Matrizen G_i)
   %A_ms = zeros(m*n, m*n);
   %for block = 0:m-2
   %   A_ms(block*n+1:(block+1)*n,     block*n+1:(block+1)*n) = reshape(G(block+1,:,:),n,n);
   %   A_ms(block*n+1:(block+1)*n, (block+1)*n+1:(block+2)*n) = - eye(n);
   %end
   %A_ms((m-1)*n+1:m*n, 1:n)           = feval(dr_y0,  Y(1,:)', Y(end,:)',param);        
   %A_ms((m-1)*n+1:m*n, (m-1)*n+1:m*n) = feval(dr_ytf, Y(1,:)', Y(end,:)',param) * G_m; 
   %---------------------------------------------
   
   % Aufbau der rechten Seite
   b_ms = zeros(n*m,1);
   b_ms(1:(m-1)*n) = - (reshape(y_ende(1:(m-1)*n), (m-1)*n, 1) - reshape(startwerte(n+1:end), (m-1)*n, 1));
   b_ms((m-1)*n+1:end) = - feval(r,Y(1,:)',Y(end,:)',param);  
   %-----------------------------------
   
   
   % Loesen des MS-Gleichungssystems zur Startwert-Korrektur
   delta_s = A_ms_sparse \ b_ms;
   
   alpha=1;
   % Neue Startwerte
   startwerte = startwerte + alpha*reshape(delta_s,n,m);

end;   


% Funktion für interne numerische Differentiation
function dz_dt = h(t, z, f, fy, n,param)
  % Vektor z zerlegen in Komponenten in Vektor y und Matrix G
  y = z(1:n);
  G = reshape(z(n+1:end),n,n);
  %G = zeros(n,n);
  %for i = 1:n
  %    G(:,i) = z(i*n+1 : (i+1)*n);
  %end;

  % DGLn
  dy_dt = feval(f, t, y,param);
  dG_dt = feval(fy, t, y,param)*G;     
  
  % Zerlegung rueckgaengig machen
  dz_dt = zeros(n*n+n,1);
  dz_dt(1:n) = dy_dt;
  dz_dt(n+1:end)=reshape(dG_dt,n*n,1);
  return;