function T12_Beispiel_Massenmatrix
%% Lösen eines DGL-Systems 2. Ordnung mit Hilfe einer Massenmatrix
%  x1''+x1'*x2'' = 1/2*(x2'+x2)
%  5*t*x1''+3*x2'' = x1-t+3*exp(t)
%  AB: x1(0)=0,x1'(0)=1,x2(0)=1,x2'(0)=1

% Anfangsbedingungen
x10   = 0; 
x10_p = 1;
x20   = 1;
x20_p = 1;
% Zeitspanne
tspan=[0 2];

%% num Integration
options=odeset('Mass', @Massenmatrix);
%                         Reihenfolge AW korrigiert (10.1.12)
[t,z]=ode45(@dgl,tspan,[x10 x20 x10_p x20_p],options);

options=odeset('Mass', @Massenmatrix, 'MaxStep',0.01);
%                         Reihenfolge AW korrigiert (10.1.12)
[ta,za]=ode45(@dgl,tspan,[x10 x20 x10_p x20_p],options);

%% Plot
subplot(2,1,1)
plot(t,z(:,1:2),'-o',t,t,t,exp(t),ta,za(:,1:2),'-o'); 
legend('x1', 'x2', 'exakte Lösung x1=t', 'exakte Lösung x2=exp(t)', ...
'x1: mit Zeitschrittbeschränkung', 'x2: mit Zeitschrittbeschränkung', ...
'Location','NorthWest'); grid on;
title('Lösungen: numerisch und exakt');
subplot(2,1,2)
plot(t,z(:,3:4),'-o',t,ones(1,length(t)),t,exp(t));
legend('x1', 'x2', 'exakte Ableitung x1´=1', 'exakte Ableitung x2´=exp(t)',...
'Location','NorthWest'); grid on;
title('Ableitungen: numerisch und exakt');
end

%% Massenmatrix
function M = Massenmatrix(t,z)
   M=[1 0  0  0 ;...
      0 1  0  0 ;...
      0 0  1  z(3); ...
      0 0 5*t 3];
end
%% Rechte Seite der DGL
function z_p=dgl(t,z)
    z_p=[z(3);... % Ergebnis: x1=t, x2=exp(t)
         z(4);...
         1/2*(z(4)+z(2));...
         z(1)-t+3*exp(t)];
end


