% c *****************************************************************
%c
%c  This program solves -u''+u=f on an interval 0<x<1 where 
%c  homogeneous Neumann boundary conditions are imposed; that is
%c  u=0 at x=0 and at x=1.  The
%c  code uses continuous piecewisequadratic  basis functions.
%c
%c *****************************************************************
%c *****************************************************************
%c
%c  INPUT
%c    nx    -  number of points in the x-direction 
%c    xl,xr -  left and right x-coordinates of domain (given interval)
%c    nnodes-  number of nodes per element or  subinterval (for linears in 1-D it's 2)
%c    nq    -  number of quadrature points per element used in the assembly
%c    nqe   -  number of quadrature points per element used in the error routine
%c
%c *****************************************************************
%c
%c  GEOMETRY
%c    xc(i)  -      x-coordinate .of global node number i 
%c    node(it,j)  - array giving global node number corresponding to
%c                  local node number j of element number it
%c    area(it)    - area of element it
%c    indx(i)     - gives unknwon number at global node number i
%c                  (i=0 if no unknown at node)
%c    xq(it,i)    - x-coordinate of quadrature point j for element it 
%c
%c    nunk        - number of unknowns
%c    nel         - number of elements
%c    np          - number of global nodes (or points)
%c *****************************************************************
%c
%c  ASSEMBLY
%c    a   -     main diagonal of the coefficient matrix 
%c    c   -     subdiagonal of coefficient matrix (matrix is symmetric)  
%c    f(i)    - right hand side array which is overwritten by solver with
%c              solution
%c
%c *****************************************************************
%c
% NOTE:  a banded Cholesky solver should be used for this problem; however
%        for simplicity we have used the built-in solver A\b
% *****************************************************************
% 
global nnodes  
global nunk  nel  neqn nq   xc  node  indx  area
global xl xr vzero kappa m czero fzero etazero 
global lambda1 lambda2 nu1 nu2  difeta beta alpha1 alpha2
global beta1 beta2 gamma1 gamma2

nx = 41;
maxstep_time = 10000;
maxstep_time = 1;
time_final = 0.1;
maxiter = 20;
maxiter = 1;
time = 0.0;
nnodes = 3;
neqn = 4;
tolerance = 0.5e-6;
%
%  Set the biological constants.
%
bio_constants 
%
%  Set up the geometry arrays.
%
geom
%
%  Assemble the mass matrix.
%
assemble_mass 
%
%  Factor the mass matrix.
%
a_mass = chol(a_mass);
% 
%  Set the initial conditions.
%
set_ic  

vcur=vold;
ccur=cold;
fcur=fold;
etacur= etaold;
%
%  TIME STEPPING LOOP
%
dx=(xr-xl)/(nx-1);
dt = 0.5 * dx * dx 
%
for  nt = 1 : maxstep_time
    
  time = time + dt
  if ( time >= time_final )
    break;
  end
%
%  ITERATION LOOP at each TIMESTEP
%
  for niter = 1 : maxiter
%
%  solve v equation
%
    disp ( 'Solve v' )
    nuk=1; 
    profile on
	vcur = solve_eqn(nuk,time,dt, a_mass, vcur,ccur,fcur,etacur,vold);  
    profile report
    profile plot
    profile off
    pause
%
%  solve c equation
%
    disp ( 'Solve c' )
    nuk=2; 
    ccur = solve_eqn (nuk,time,dt, a_mass, vcur,ccur,fcur,etacur,cold);   
%
%  solve f equation
%
    disp ( 'Solve f' )
    nuk=3; 
	fcur = solve_eqn (nuk,time,dt,a_mass, vcur,ccur,fcur,etacur,fold);  
%
%  solve eta equation
%
    disp ( 'Assemble stiff' )
    nuk=4; 
	assemble_stiff  
    disp ( 'Factor stiff' )
	a_stiff = chol(a_stiff);
    disp ( 'Solve eta' )
	etacur = solve_eqn (nuk,time,dt,a_stiff, vcur,ccur,fcur,etacur,etaold);   
%
% calculate residual; if less than tolerance, go to next time step
%
    disp ( 'Compute residual' )
    residual   
    if ( resid_norm <= tolerance)
      break;
    end 
%
% if failed to converge in maxiter steps, then reduce time step and start again
%
    if ( niter == maxiter )
      disp ( 'Reducing time step.' );
	  time = time-dt;
	  dt = dt/2;
	  vcur=vold;
      ccur=cold;
      fcur=fold;
      etacur= etaold;
	end
    

  end
%
% get ready for next time step
%
  vold = vcur;  
  cold = ccur;
  fold = fcur;
  etaold = etacur;

end