function out = ex12ode(t,y,flag)

% Copyright 1999, The MathWorks, Inc.

global P N

if nargin < 3 | isempty(flag)
  out = f(t,y,P,N);
else
  switch(flag)           
  case 'mass'
    %  The first N equations approximate the elliptic PDE so are algebraic
    %  in nature.  The next N equations approximate the parabolic PDE so are
    %  ODEs except for the first which is algebraic because it comes from a
    %  Dirichlet boundary condition.
    out = diag([zeros(N+1,1); ones(N-1,1)]);
  case 'jpattern'
    out = Jacpat(N);
  otherwise
    error(['Unknown flag ''' flag '''.']);
  end
end

%--------------------------------------------------------------
function yp = f(t,y,P,N)
%  Evaluate the differential equations.

u1 = y(1:N);
u2 = y(1+N:N+N);
f1 = zeros(N,1);
f2 = zeros(N,1);

%  The are N equally spaced mesh points in 0 <= x <= 1:
h = 1/(N - 1);

%  Equations at interior mesh points:
for i = 2:(N-1)
  f1(i) = (u1(i+1) - 2*u1(i) + u1(i-1))/h^2 + exp(-2*u2(i));              
  f2(i) = (u2(i+1) - 2*u2(i) + u2(i-1))/h^2 + exp(-2*u2(i)) + exp(-u1(i));
end

%  Boundary conditions at x = 0:
f1(1) = ((u1(2) - u1(1))/h - exp(-u1(1)))/h - exp(-2*u2(1));
f2(1) = u2(1) - log(P + t);

%  Boundary conditions at x = 1:  
f1(N) = u1(N) - log(1 + P + t);
f2(N) = (exp(-u1(N)) - (u2(N) - u2(N-1))/h)/h + exp(-2*u2(N)) + exp(-u1(N));

yp = [f1; f2];

%--------------------------------------------------------------
function S = Jacpat(N)
%  Define the sparsity structure of the Jacobian by indicating which 
%  variables appear in which equations.  This is done by inspection of
%  the function 'f' above for evaluating the equations.

S = sparse(2*N,2*N);
S(1,1) = 1;          %  u1(1)
S(1,2) = 1;          %  u1(2)
S(1,1+N) = 1;        %  u2(1)
for i = 2:(N-1)
  S(i,i-1) = 1;      %  u1(i-1)
  S(i,i) = 1;        %  u1(i)
  S(i,i+1) = 1;      %  u1(i+1) 
  S(i,i+N) = 1;      %  u2(i)
end
S(N,N) = 1;          %  u1(N)  
S(1+N,1+N) = 1;      %  u2(1)
for i = 2:(N-1)
  S(i+N,i) = 1;      %  u1(i)
  S(i+N,i-1+N) = 1;  %  u2(i-1)
  S(i+N,i+N) = 1;    %  u2(i)
  S(i+N,i+1+N) = 1;  %  u2(i+1)
end
S(N+N,N) = 1;        %  u1(N)
S(N+N,N-1+N) = 1;    %  u2(N-1)
S(N+N,N+N) = 1;      %  u2(N)  
