function Ufinal = vera(J,N,tf,Uinit)
% VERA Solves the PDE   u_t = (1+2x) u_xx  better than BOB does.
% Result is plotted and final approximation of  u(x,t_f)  is returned.
% Example of usage with default u(x,0) = x(1-x):
%   >> J=20; N=1000; tf=0.1; UN = vera(J,N,tf);
% Examples showing behavior at critical value of mu (for stability):
%   >> J=100; N=5690; tf=0.1; UN = vera(J,N,tf);
%   >> J=100; N=5685; tf=0.1; UN = vera(J,N,tf);
% Example of usage with initial u(x,0) given as fourth argument:
%   >> J=25; x=0:1/J:1;  U0=zeros(size(x));
%   >> U0(1:J/5)=1;  U0=U0+fliplr(U0);  U0
%   >> vera(J,200,0.03,U0);
% ELB 2/4/07

dx = 1/J;  mu = (tf/N) / (dx*dx);  % = dt/dx^2
disp(sprintf('vera is using mu = %f',mu))
x = 0:dx:1;
if nargin < 4,  Uinit = x .* (1 - x);  end

U = Uinit;
for n=1:N
   U(1)=0; U(J+1)=0;
   U(2:J) = U(2:J) + mu * (1+2*x(2:J)) .* (U(3:J+1) - 2 * U(2:J) + U(1:J-1));
   if any( (U==NaN) | (U==inf) | (U==-inf) )
      error('vera:blowup','Solution has exploded by step %d.  Stopping.',n)
   end
end
Ufinal = U;

if J < 50, plot(x,Uinit,x,Ufinal,'o--'); 
else, plot(x,Uinit,x,Ufinal,'--'); end
legend('u(x,0)','U_j^N ~~ u(x,t_f)'), xlabel x, ylabel u