function [sol, I] = transistormol(n)
%
% transistor equation from pdex3 by the method of lines
%
close all
% Problem parameters, shared with nested functions.
L = 1;
d = 0.1;
eta = 10;
K = d;
I_p = 1;
% pdex3 has 0 <= x3 <= 1.  Here we use -1 <= x <= 1.
% Change of variable means x = 2x3 - 1, hence d/dx = (1/2) d/dx3. 
D = gallery('chebspec',n+2);
x = cos( pi*(0:n+1)/(n+1) );
function dy = mol(t,y)
    my = [0; y; 0];
    % tmp = -eta/L*2*D*my + d*4*D*D*my  is the comment version pde
    tmp = -d*eta*2*D*my + d*4*D*D*my; % is the coded version pdex3
    dy = [tmp(2:end-1)];
end
x3 = (x+1)/2;
% (K*L/d)*(1 - exp(-eta*(1 - x)))/eta
inits = K*L/d/eta*(1-exp(-eta*(1-x3/L)));
tf = 1.0;
sol = ode15s( @mol, [0,tf], inits(2:end-1)' );
t = linspace(0,tf,101);
z = [zeros(size(sol.x)); sol.y ; zeros(size(sol.x)) ];
zx = 2*D*z;
figure
I = I_p*d*zx(end,:)/K;
semilogy( sol.x, I, 'k.')
xlabel('t')
ylabel('Emitter discharge current I(t)')
figure
%plot( x, inits, 'k+', x , sol.y(:,end), 'ko' )
mesh( x3(2:end-1), sol.x, sol.y' ), view(160,50),colormap([0,0,0]), axis([0,1,0,tf,0,0.1]);
end

