function prob_2_36
% 
% Written by:
% -- 
% John L. Weatherwax                2005-05-07
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----
  
global J h
N = 50;
h = 1/(N+1);
x = h*(1:N);

% Initial conditions: 
v0 = zeros(size(x));
ids = find( x <= 0.1 ); v0(ids) = 0;
ids = find( 0.1 < x & x < 0.25 ); v0(ids) = (x(ids)-0.1)/0.15;
ids = find( 0.25 < x & x <= 1.0 ); v0(ids) = 1.0;
e = ones(N,1);
J = spdiags([e -2*e e],-1:1,N,N)/h^2;
J(end,end-1) = 2*J(end,end-1);
options = odeset('Jacobian',J);
tspan = 0:0.001:0.006;
%tspan = 0:0.1:0.6;
tic; [t,v] = ode15s(@f,tspan,v0,options); toc
tic; [t,v] = ode23(@f,tspan,v0,options); toc
% Add the boundary values:
x = [0 x 1];
npts = length(t);
v = [zeros(npts,1) v v(:,end-1)];
mesh(x,t,v)
xlabel( 'x' ); ylabel( 't' ); 

%=========================================
function dvdt = f(t,v)
global J h
dvdt = J*v + g(t,v);
%dvdt = J*v;

%=========================================
function res = g(t,v)
  
inds1 = find( v <= 0.5 ); 
res(inds1,:) = -(2*10^5)*v(inds1,:);

inds2 = find( v > 0.5 );   
res(inds2,:) = 0;