function ex8
%  This is the two transistor amplifier problem of E. Hairer, C. Lubich, 
%  and M. Roche, The Numerical Solution of Differential-Algebraic Systems
%  by Runge-Kutta Methods, Lecture Notes in Math. #1409, Springer, Berlin,
%  1989.  It is solved in the form M*u' = f(t,u) with a mass matrix that is
%  not diagonal.  The mass matrix is treated as a sparse matrix to test the
%  capability.  Fig. 9.4 of HLR shows the solution on [0 0.2], but here it
%  is computed on [0 0.05] because the computation is less expensive and the
%  nature of the solution is clearly visible on the shorter interval.

% Copyright 2004, The MathWorks, Inc.

%  Hairer and Wanner's RADAU5 requires consistent initial conditions
%  which they take to be
Ub = 6;
u0 = zeros(8,1);
u0(2) = Ub/2;
u0(3) = u0(2);
u0(4) = Ub;
u0(5) = Ub/2;
u0(6) = u0(5);
u0(7) = Ub;
u0(8) = 0;

% Create Mass Matrix
c = 1e-6 * [1:5];
M = sparse(8,8);
M(1,1) = -c(1);
M(1,2) =  c(1);
M(2,1) =  c(1);
M(2,2) = -c(1);
M(3,3) = -c(2);
M(4,4) = -c(3);
M(4,5) =  c(3);
M(5,4) =  c(3);
M(5,5) = -c(3);
M(6,6) = -c(4);
M(7,7) = -c(5);
M(7,8) =  c(5);
M(8,7) =  c(5);
M(8,8) = -c(5);
    
options = odeset('Mass',M,'MassSingular','yes');

[t,u] = ode15s(@ex8ode,[0 0.05],u0,options);

Ue = 0.1*sin(200*pi*t);

figure(1)
plot(t,Ue,t,u(:,8))
legend('input ','output ',3)
title('Two Transistor Amplifier Problem')
axis([0 0.05 -5 1])

%-------------------------------------------------------------------------%

function yp = ex8ode(t,u)

  Ub = 6;
  R0 = 1000;
  R = 9000;
  alpha = 0.99;
  beta = 1e-6;
  Uf = 0.026;

  Ue = 0.1*sin(200*pi*t);
  f23 = beta*(exp((u(2) - u(3))/Uf) - 1);
  f56 = beta*(exp((u(5) - u(6))/Uf) - 1);
  f = zeros(8,1);
  f(1) = - ((Ue - u(1))/R0);
  f(2) = - ((Ub - u(2))/R - u(2)/R + (alpha - 1)*f23);
  f(3) = - (f23 - u(3)/R);
  f(4) = - ((Ub - u(4))/R - alpha*f23);
  f(5) = - ((Ub - u(5))/R - u(5)/R + (alpha - 1)*f56);
  f(6) = - (f56 - u(6)/R);
  f(7) = - ((Ub - u(7))/R - alpha*f56);
  f(8) = - (- u(8)/R);
  yp = f;
