function ex7
%  This is the one transistor amplifier problem displayed on p. 377 of 
%  E. Hairer and G. Wanner, Solving Ordinary Differential Equations II
%  Stiff and Differential-Algebraic Problems, 2nd ed., Springer, Berlin,
%  1996.  It is easy enough to write this problem in semi-explicit form,
%  but it arises in the form M*u' = f(t,u) with a mass matrix that is
%  not diagonal and it is solved here in the form it arises.  The mass
%  matrix is treated as a sparse matrix to test the capability. Fig. 1.4
%  of Hairer and Wanner 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 2002, The MathWorks, Inc.

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

%  Perturb the algebraic variables to test initialization.  Note that
%  the initial values computed here approximate the ones used by Hairer
%  and Wanner, but they are different.
u0(4) = u0(4) + 0.1;
u0(5) = u0(5) + 0.1; 

%create Mass Matrix
c = 1e-6 * [1:3];
M = sparse(5,5);
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);
options = odeset('Mass',M);

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

Ue = 0.4*sin(200*pi*t);
figure(1)
plot(t,Ue,t,u(:,5))
legend('input ','output ',2)
title('One Transistor Amplifier Problem')

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

function yp = ex7ode(t,u)

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

  Ue = 0.4*sin(200*pi*t);
  f23 = beta*(exp((u(2) - u(3))/Uf) - 1);
  f = zeros(5,1);
  f(1) = - (Ue - u(1))/R0;
  f(2) = - (Ub/R - u(2)*2/R - 0.01*f23);
  f(3) = - (f23 - u(3)/R);
  f(4) = - ((Ub - u(4))/R - 0.99*f23);
  f(5) = - (- u(5)/R);
  yp = f;
