Re: Colnew Help

/* This is Example 4 with the notations as in  http://arxiv.org/pdf/hep-th/0407005 */

y:1.9d0$ aleft : 0.0d0$ aright : y$ ntol:4$
depends([p,f2,f3],t)$
lvfs:[t,p,f2,f3]$
list_of_eqns: [diff(p,t,2)= -0.5*(1/t+1/(t-1.0)+1/(t-y))*diff(p,t,1)+
((v-f2-f3)/t+f2/(t-1.0)+f3/(t-y))*p , diff(f2,t,1)=0, diff(f3,t,1)=0]$
list_of_conds:[(diff(p,t,1)-2.0*p*(v-f2-f3))(0)=0,
(p -1.0)(1)=0, (diff(p,t,1)-2.0*p*f2)(1)=0, (diff(p,t,1)-2.0*p*f3)(y)=0]$
list_of_guess:[2.0*t-1.0,1.0,1/(2.0*y-1.0)]$


load(colnew)$
load("colnew_discov.lisp")$

print("Define the interval [aleft,aright] for the differential equations.");

if (aleft = 'aleft) then
aleft: float(read("aleft: ")); /* terminate answer with ; */
if (aright = 'aright) then
aright: float(read("aright: "));

print("Fix the number ntol of checked tolerances");
if (ntol = 'ntol) then ntol: read("ntol: "); 

print("Give the list of variable and functions in equations, like [x,u,v]");

if (lvfs='lvfs) then lvfs :read("lvfs = ");
depends(rest(lvfs),first(lvfs))$

print("Define the list of the differential equations.");
  /* Each one of the form
  u^(n)=f(lower order terms), where depends([u,...],var0). Example:
  if {u''= u'+u*v,v'=-v+2*u'} then depends([u,v],x)$
  leqs:[diff(u,x,2)=diff(u,x,1)+u*v,diff(v,x,1)=-v+2*diff(u,x,1)]$
  The equations may depend on extra parameters
  which may be varied for continuation, see Example 3. */

if (list_of_eqns = 'list_of_eqns) then
list_of_eqns: read("list_of_eqns: "); /* list-of-eqns" */
eqns_to_funs(list_of_eqns)$

/* Now max-subst-list mstar ncomp m-list var-list var0 and f are defined,
   and expressed on the z-vector, that is the list of mstar lower order
   derivatives. */

define(funmake(df,cons(var0,var_list)), jacobian(apply(f,cons(var0,var_list)),var_list))$

print("Define the ", mstar, "boundary conditions.");
  
  /* Each one of the form C(u, u',...)(p)=0 where only lower order
     derivatives appear and the condition is at point p, which may be aleft, aright,
     or an interior point. */

  /* The interior points form an array, ordered, named fixpnt, of size nfixpnt.
     Example: (u-5)(0)=0, (u'+2)(1)=0, (v-1)(0.5)=0. Note 0.5 above, 1/2 does not work. */

if (list_of_conds = 'list_of_conds) then
list_of_conds: read("list_of_conds: ");  /* list-of-conds  boundary conditions */
conds_to_zconds(list_of_conds, aleft, aright, max_subst_list, m_list)$
fixpnt: float(fixpnt)$
ZETA:float(ZETA)$  


  /* Now g, ZETA, fixpnt are defined. ZETA is the list of length mstar of
     points at which boundary conditions are defined, fixpnt the list of the
     interior ones. One may add further points to fixpnt, for example to avoid
     singularities. Function g is expressed on the z-vector. */

define(funmake(dg, var_list), jacobian(apply(g,var_list),var_list))$


   /* In general the considered systems have several discrete solutions
      fulfilling both equations and boundary conditions. To select and converge
      on a particular one, it is necessary to provide a guess as starting point.
      This is a list of functions of var0 corresponding to the ncomp underived
      functions of the system. The following uses maxima to compute their
      appropriate derivatives. */

if (list_of_guess = 'list_of_guess) then
list_of_guess: read("list_of_guess = ");
define(guess(var0), build_guess(list_of_guess,m_list,ncomp,var0))$


   /* Finally one must build the list ipar which controls the running
      of colnew.  ipar[1] is set to one for nonlinear problem. In the much
      simpler linear case set ipar[1] to 0. One assumes there is a guess so
      we set ipar[9] to 1. If we use continuation, that is one uses a
      previous solution as guess to a new one after having varied some
      parameters, one needs to set manually ipar[9]=3 and
      ipar[3]=ispace[1]. See Example 3. Also by default one sets tolerances
      1.0D-5 on all components of the z-vector, if some are to be excluded,
      one needs to change by hand ntol, tol, and ltol. In this case one can
      assign ntol: <> before then tol will have ntol length, and will have
      to be edited. */

if ntol='ntol then ntol: mstar;
build_ipar(m_list, mstar, ncomp, nfixpnt,ntol)$
fspace:makelist(0d0,k,1,ipar[5])$
ispace:makelist(0,k,1,ipar[6])$


   /* Better check the colnew variables evaluated before. One can still modify
      them, notably tol and ltol or fixpnt. */

print("ncomp = ", ncomp, ", m_list = ", m_list, ", mstar = ", mstar,
", ntol = ", ntol, ", max_subst_list = ", max_subst_list, "");
print("ZETA: ", ZETA);
print("ipar: ", ipar);
print("ltol: ", ltol);
print("tol: ",tol);
print("fixpnt: ",fixpnt);
print("nfixpnt: ",nfixpnt);


/* At this point all parameters in f, df, g, dg must be fixed to numerical values.
   In general there may still contain free parameters used for continuation. One must
   give here values to these parameters, and define fsub, dfsub, gsub, dgsub whose purpose
   is to evaluate f, df, g, dg with these parameters. Moreover one should start with
   iflag=1 or 0 in the singular case. In the case of continuation, one has to run
   the following in a loop, see Example 3, while varying parameters and modifying ipar[9].
   In particular iflag, fspace, ispace will be modified in the loop. */

iflag:1$
print("Please give floating numerical values to parameters.");
/* ipar[7]:1$   to suppress abundant colnew printout */

ipar[7]:1$

break();

/* Example3: Print ltol:[1,3]$ tol:[1d-7,1d-7]$ s:0.2d0$ L:60d0$ followed by exit; */

nsteps:100$

/*  We store f2, f3 and the energy f2+yf3 in the array eigenval. Since these are constants,
we compute them with colnew_appsln at t=0.1d0 */

eigenval: makelist([0.0,0.0,0.0],k,1,nsteps)$

/* We run the continuation computation in a loop. */

v:0.00001d0$
for ii from 1 thru nsteps do (
define(funmake (fsub, cons(var0,var_list)), apply (f, cons(var0,var_list))),
define(funmake (dfsub, cons(var0,var_list)), apply(df,cons(var0,var_list))),
define(funmake (gsub, cons(i,var_list)), buildq([gval:apply (g, var_list)], gval[i])),
define(funmake (dgsub, cons(i,var_list)), buildq([dgval:apply(dg,var_list)], row(dgval,i)[1])),
[iflag, fspace, ispace] : colnew_expert(ncomp, m_list, aleft, aright, ZETA,
  ipar, ltol, tol, fixpnt, ispace, fspace, iflag, fsub, dfsub, gsub,
  dgsub, guess),
zval: colnew_appsln([0.1d0], 4, fspace, ispace)[1],
eigenval[ii][1]: zval[3],  /* zval[1] zval[2] are p and p', useless. */
eigenval[ii][2]: zval[4],
eigenval[ii][3]: zval[3]+y*zval[4],
ipar[9]:3,
ipar[3]: ispace[1],
v:v+0.2)$ /* set continuation after first step */

time(%);

vv: makelist(0.2*k,k,0,nsteps-1)$
plot2d([[discrete,vv, makelist(eigenval[k][1],k,1,nsteps)],
  [discrete,vv, makelist(eigenval[k][2],k,1,nsteps)],
  [discrete,vv, makelist(eigenval[k][3],k,1,nsteps)]])$

;;; Input in the form:  diff(F,x,n)=... (lower order terms), conditions in the form
;;; C(F,...)(point)=0                     . Example (see example3):
;;; depends([G,H],x)$ 
;;; list-of-eqns:[ diff(G,x,2)=L^2*s*(G-1)-L*((3-n)/2*H*diff(G,x)+(n-1)*diff(H,x)*G),
;;; diff(H,x,3)= L^3*(1-G^2)+L^2*s*diff(H,x)-L*((3-n)/2*H*diff(H,x,2)+n*(diff(H,x))^2)]$
;;; list-of-conds:[G(0),H(0),diff(H,x)(0),(G-1)(1),diff(H,x)(1)]$
;;; Further implicit parameters, like L, s, n above may appear and used in continuation
;;; Output: n, variables z1,...,zn, and function f(x,z1,....,zn), a list of all eqs.
;;; Idem for conditions a list g(z1,....,zn) of all conds on cond points.

;;; The functions f(x,z1,....,zn) and g(z1,....,zn) are stored as maxima functions, hence
;;; define special variables $f, $g to which they are attached.
;;; Also give access to the many parameters of colnew: ncomp, mstar, ntol etc.

(declaim (special $f $g  $zvec $dmval $var0 $var_list $m_list $ntol $nfixpnt $fixpnt
		  $ncomp $mstar $aleft $aright $ipar |$zeta|))

(defun eq-disc (eqn)

  "Treat one maxima equation. It should be of the form diff(f,x,n)=... Check that."

  (if (not (equal (caar eqn) 'MEQUAL))
      (merror (intl:gettext "eqn:  needs an equation.")))
  (let* ((lhs (cadr eqn))
	 (diffop (caar lhs))
	 (functn (cadr lhs))
	 (variable (caddr lhs))
	 (order (cadddr lhs)))
    (if (not (equal diffop '%DERIVATIVE))
	(merror (intl:gettext "eqn: needs a differential equation.")))
    (if (> order 4) (merror (intl:gettext "Maximum degree of equation 4")))
    (values order functn variable)))

(defun cond-disc (condt)
  
  "Treat one cond, extract the condition and the point at which it is given. To keep things simple
   assume the condition is given in the form C(F,G,..)(p)=0 at point p. For example:
   :lisp #$$ (G+diff(G,x)-1)(3)=0$ -> ((MEQUAL) ((MQAPPLY) ((MPLUS) $g (($DIFF) $g $X) ((MMINUS) 1)) 3) 0)"
  
  (if (not (and (eql (caar condt) 'MEQUAL) (eql (caddr condt) 0)))
      (merror (intl:gettext "condt: must be of the form C(p)=0.")))
  (let* ((lhs (cadr condt)) 
	 (condit (if (eql (caar lhs) 'MQAPPLY) (cadr lhs) (caar lhs)))
	 (point (if (eql (caar lhs) 'MQAPPLY) (caddr lhs) (cadr lhs))))
    (values condit point)))


    
(defun $eqns_to_funs (list-of-eqns)

  "Treat a maxima list of maxima equations and a list of constraints. Outputs ncomp, mstar,
    m-list, for colnew, and assigns maxima function f and several maxima variables."

  (if (not (equal (caar list-of-eqns) 'MLIST))
      (merror (intl:gettext "Needs a list of equation.")))
  (let ((mstar 0) (ncomp 0) (index 0) m-list  newvar subst-list
	var0 var-list funs-list)
    (setq list-of-eqns (cdr list-of-eqns)) ; from maxima list to lisp list

;;; Treat equations
    (dolist (eqn list-of-eqns)
      (incf ncomp)
      (multiple-value-bind (order functn variable) (eq-disc eqn)
	(if (= 0 index) (setf var0 variable)
	    (if (not (equal variable var0))
		(merror (intl:gettext "Only one variable!"))))
	(push order m-list)
;;;  Introduce the Z-vector
	(dotimes (i order)
	  (setq newvar (intern (concatenate 'string "$Z" (write-to-string (incf index)))))
	  (push newvar var-list)
	  (if (= i 0)
	      (push `((MEQUAL)  ,functn ,newvar) subst-list)
	      (push `((MEQUAL)
		      ((%DERIVATIVE) ,functn ,variable ,i) ,newvar)
		    subst-list)))))
    (setf var-list (nreverse var-list))
    (if (> ncomp 20) (merror (intl:gettext "Maximum number of equations 20")))
    (setf mstar (reduce #'+ m-list))
    (setq m-list (cons '(MLIST) (nreverse m-list)))
    (if (> mstar 40) (merror (intl:gettext "Maximum sum of degrees 40")))
    (let ((max-subst-list (cons '(MLIST) subst-list)))
      (dolist (eqn list-of-eqns)
	; Doesn't work if max-subst-list is reversed!
	(push (simplify ($substitute max-subst-list (caddr eqn))) funs-list))
      (setq funs-list (nreverse funs-list))    
;;; Export variables to maxima
      (meval  `((MSETQ) $max_subst_list ,max-subst-list)))
    (meval  `((MSETQ) $mstar ,mstar))
    (meval  `((MSETQ) $ncomp ,ncomp))
    (meval  `((MSETQ) $m_list ,m-list))
    (meval  `((MSETQ) $var_list ,(cons '(MLIST) var-list)))
    (meval  `((MSETQ) $var0 ,var0))
    (meval  `((MDEFINE) ,(cons '($f) (list* var0  var-list)) ,(cons '(MLIST) funs-list))))
  'DONE)

(defun $conds_to_zconds(list-of-conds $aleft $aright $max_subst_list $m_list)
  
  "Express conditions in terms of the zvector, put them in order and discover if there are interior points to be
   inserted into the list of fix points."
  
  (setq list-of-conds (cdr list-of-conds))
  (if (not (eql (length list-of-conds) $mstar))
      (merror (intl:gettext "Need to have as many conditions as free parameters")))
  (let ((conds-list ()) (fixpnt ()) (zet ()) nfixpnt)
    (dolist (condt list-of-conds)
      (multiple-value-bind (condit point) (cond-disc condt)
	(push (list point (simplify ($substitute $max_subst_list  condit))) conds-list)))
    (setq conds-list (nreverse conds-list))
;;; Ensure points of ZETA are in order and constraints are sorted accordingly
    (setq conds-list (stable-sort conds-list #'< :key #'car))
    (setq zet (mapcar #'car conds-list))
    (setq conds-list (mapcar #'cadr conds-list))
    (dolist (zz zet) (if (and (> zz $aleft) (< zz $aright)
			      (not (member zz fixpnt)))
			 (push zz fixpnt)))
    (setq fixpnt (nreverse fixpnt))
    (setq nfixpnt (length fixpnt))
;;; Define  maxima function g, and maxima variables fixpnt ZETA (zeta is the function zeta)
    (meval  `((MDEFINE) ,(cons '($g) (cdr $var_list)) ,(cons '(MLIST) conds-list)))
;;; |$zeta| becomes ZETA in maxima (Manual, section 37.1) Other maxima variables
    (meval  `((MSETQ) |$zeta| ,(cons '(MLIST) zet)))
    (meval  `((MSETQ) $fixpnt ,(cons '(MLIST) fixpnt)))
    (meval  `((MSETQ) $nfixpnt ,nfixpnt)))
  'DONE)
    


(defun $build_guess ($list_of_guess $m_list $ncomp $var0)

  "list-of-guess is the init-guess, a maxima list of ncomp functions, giving a guess for the
   underivated functions. Construct z-vector and dmval in terms of var0, the variable of diff. eqns."

  (if (not (and (eql (caar $list_of_guess) 'MLIST)
		(eql ($length $list_of_guess) $ncomp)))
      (merror (intl:gettext "Guess must be a maxima list of functions.
                             of length ncomp")))
  (let (z-vec  dmval (ml (cdr $m_list)))
    (dolist (gfun  (cdr $list_of_guess))
      (let ((gexp gfun))
;;; For each guess compute derivatives of guess as indicated by m-list
	(dotimes (k (car ml))
	  (push gexp z-vec)
	  (setf gexp (simplify (meval `(($diff) ,gexp ,$var0)))))
	(push gexp dmval)
	(setq ml (cdr ml))))
    (setq z-vec (nreverse z-vec))
    (setq dmval (nreverse dmval))
    (setq $zvec (cons '(MLIST) z-vec))
    (setq $dmval (cons '(MLIST) dmval))
    (list '(MLIST) $zvec $dmval)))


;;; Example:  list_of_guess: [cos(x) + x^2,x^5]$ m_list:[2,3]$ ncomp:2$ var0:x$
;;; result: [[cos(x)+x^2,2*x-sin(x),x^5,5*x^4,20*x^3],[2-cos(x),60*x^2]]


(defun $build_ipar ($m_list $mstar $ncomp $nfixpnt $ntol &optional ($nmax 40) ($nsubint 5))
  
  ;; by default nmax=40 ntol=mstar tolerances set to 1d-5 ncomp=5
  ;; ntol is the number of tolerances, ltol specifies to what Z components they apply, tol their values.
  ;; for Problem 3 ipar=[1,4,10,2,40000,2500,1,0,1,0,0]
  ;; nmax is the maximum number of subintervals, ispace and fspace are resp  nmax*nsizei  nmax*nsizef.
  ;; One can add further fixed points in the array fixpnt and increase accordingly nfixpnt for example to
  ;; avoid a singularity in the equations, since the equations are never applied on fixed points. See Problem 4.
  
  (let ((mmax (apply 'max (cdr $m_list))))
    (if (= $ncomp 0)  (setq $ncomp (max (1+ mmax) (- 5 mmax)))
	(progn (if (< $ncomp mmax) (setq $ncomp mmax))
	       (if (> $ncomp 7) (setq $ncomp 7))))
    (if (= $nsubint 0) (setq $nsubint 5))
    (if (< $ntol 1) (setq $ntol 1))
    (if (> $ntol $mstar) (setq $ntol $mstar))
    (prepare-tol $ntol $mstar) ; fill tol and ltol with basic defaults
    (let* ((kd (* $nsubint $ncomp))
	   (kdm (+ kd $mstar))
	   (nsizei (+ 3 kdm))
	   (nsizef (+ 4 (* 3 $mstar) (* kdm (+ 5 kd)) (* 4 $mstar $mstar)))
	   ;; space is cheap one sets  (2*mstar-nrec)-> 2*mstar and also finally
	   (ispace (max (* $nmax nsizei) 2500))
	   (fspace (max (* $nmax nsizef) 40000))
	   ;; Here one sets ipar[9]=1 assuming there is a guess, for continuation one has to set
	   ;; ipar[9]=3 and ipar[3]:ispace[1] see the example in Problem 3.
	   (ipar (list '(MLIST) 1 $ncomp $nsubint $ntol fspace ispace 0 0 1 0 $nfixpnt)))

      (setq $ipar ipar)))
  'DONE)		


(defun prepare-tol ($ntol $mstar)
  
  ;; ltol is filled with increasing indices k (1<=k<=mstar) for which the tolerance has to be checked.
  ;; tol gives the corresponding tolerances. For example in problem 3 one has mstar=5, but ntol=2, and one chooses
  ;; ltol=[1,3], so Z1 and Z3 are checked, and tol=[1d-7,1d-7]. To see the correspondance with original variables
  ;; print max_subst_list, here only G and H are checked. If the tolerance on Z_k is eps, that means that
  ;; |Z_k - Z_k^0| < eps * Z_k^0 + eps where Z_k^0 is the exact solution and of course this difference is estimated
  ;; by looking at the rate of convergence of the Z_k after successive subdivisions.
  
  (defparameter $ltol  (cons '(MLIST) (loop for j below $ntol collect (1+ j))))
  (defparameter $tol (cons '(MLIST) (loop for f below $ntol collect 1.0D-5)))
  (if (< $ntol $mstar) (format t "You should edit the arrays ltol and tol ~%")
      (format t "The array tol is filled by default, you can edit tol ~%"))
  'DONE)

load(colnew)$
load("colnew_discov.lisp")$

/* Data corresponding to Example 3, first iteration. */

/* n:0.2$ cons:0.5*(3.0-n)$ ntol:2$ aleft:0.0$ aright: 1.0$
lvfs: [x,G,H]$
depends([G,H],x)$
list_of_eqns:
[ diff(G,x,2)=L^2*s*(G-1)-L*((3-n)/2*H*diff(G,x)+(n-1)*diff(H,x)*G),
diff(H,x,3)= L^3*(1-G^2)+L^2*s*diff(H,x)-L*((3-n)/2*H*diff(H,x,2)+n*(diff(H,x))^2)]$
list_of_conds:
[G(0)=0,(G-1)(1)=0,H(0)=0,diff(H,x)(1)=0,diff(H,x)(0)=0]$
list_of_guess:[1.0-exp(-L*x),-L^2*x^2*exp(-L*x)]$ */

/* End of data */

print("Define the interval [aleft,aright] for the differential equations.");

if (aleft = 'aleft) then
aleft: float(read("aleft: ")); /* terminate answer with ; */
if (aright = 'aright) then
aright: float(read("aright: "));

print("Fix the number ntol of checked tolerances");
if (ntol = 'ntol) then ntol: read("ntol: "); 

print("Give the list of variable and functions in equations, like [x,u,v]");

if (lvfs='lvfs) then lvfs :read("lvfs = ");
depends(rest(lvfs),first(lvfs))$

print("Define the list of the differential equations.");
  /* Each one of the form
  u^(n)=f(lower order terms), where depends([u,...],var0). Example:
  if {u''= u'+u*v,v'=-v+2*u'} then depends([u,v],x)$
  leqs:[diff(u,x,2)=diff(u,x,1)+u*v,diff(v,x,1)=-v+2*diff(u,x,1)]$
  The equations may depend on extra parameters
  which may be varied for continuation, see Example 3. */

if (list_of_eqns = 'list_of_eqns) then
list_of_eqns: read("list_of_eqns: "); /* list-of-eqns" */
eqns_to_funs(list_of_eqns)$

/* Now max-subst-list mstar ncomp m-list var-list var0 and f are defined,
   and expressed on the z-vector, that is the list of mstar lower order
   derivatives. */

define(funmake(df,cons(var0,var_list)), jacobian(apply(f,cons(var0,var_list)),var_list))$

print("Define the ", mstar, "boundary conditions.");
  
  /* Each one of the form C(u, u',...)(p)=0 where only lower order
     derivatives appear and the condition is at point p, which may be aleft, aright,
     or an interior point. */

  /* The interior points form an array, ordered, named fixpnt, of size nfixpnt.
     Example: (u-5)(0)=0, (u'+2)(1)=0, (v-1)(0.5)=0. Note 0.5 above, 1/2 does not work. */

if (list_of_conds = 'list_of_conds) then
list_of_conds: read("list_of_conds: ");  /* list-of-conds  boundary conditions */
conds_to_zconds(list_of_conds, aleft, aright, max_subst_list, m_list)$
fixpnt: float(fixpnt)$
ZETA:float(ZETA)$  


  /* Now g, ZETA, fixpnt are defined. ZETA is the list of length mstar of
     points at which boundary conditions are defined, fixpnt the list of the
     interior ones. One may add further points to fixpnt, for example to avoid
     singularities. Function g is expressed on the z-vector. */

define(funmake(dg, var_list), jacobian(apply(g,var_list),var_list))$


   /* In general the considered systems have several discrete solutions
      fulfilling both equations and boundary conditions. To select and converge
      on a particular one, it is necessary to provide a guess as starting point.
      This is a list of functions of var0 corresponding to the ncomp underived
      functions of the system. The following uses maxima to compute their
      appropriate derivatives. */

if (list_of_guess = 'list_of_guess) then
list_of_guess: read("list_of_guess = ");
define(guess(x), build_guess(list_of_guess,m_list,ncomp,var0))$


   /* Finally one must build the list ipar which controls the running
      of colnew.  ipar[1] is set to one for nonlinear problem. In the much
      simpler linear case set ipar[1] to 0. One assumes there is a guess so
      we set ipar[9] to 1. If we use continuation, that is one uses a
      previous solution as guess to a new one after having varied some
      parameters, one needs to set manually ipar[9]=3 and
      ipar[3]=ispace[1]. See Example 3. Also by default one sets tolerances
      1.0D-5 on all components of the z-vector, if some are to be excluded,
      one needs to change by hand ntol, tol, and ltol. In this case one can
      assign ntol: <> before then tol will have ntol length, and will have
      to be edited. */

if ntol='ntol then ntol: mstar;
build_ipar(m_list, mstar, ncomp, nfixpnt,ntol)$
fspace:makelist(0d0,k,1,ipar[5])$
ispace:makelist(0,k,1,ipar[6])$


   /* Better check the colnew variables evaluated before. One can still modify
      them, notably tol and ltol or fixpnt. */

print("ncomp = ", ncomp, ", m_list = ", m_list, ", mstar = ", mstar,
", ntol = ", ntol, ", max_subst_list = ", max_subst_list, "");
print("ZETA: ", ZETA);
print("ipar: ", ipar);
print("ltol: ", ltol);
print("tol: ",tol);
print("fixpnt: ",fixpnt);



/* At this point all parameters in f, df, g, dg must be fixed to numerical values.
   In general there may still contain free parameters used for continuation. One must
   give here values to these parameters, and define fsub, dfsub, gsub, dgsub whose purpose
   is to evaluate f, df, g, dg with these parameters. Moreover one should start with
   iflag=1 or 0 in the singular case. In the case of continuation, one has to run
   the following in a loop, see Example 3, while varying parameters and modifying ipar[9].
   In particular iflag, fspace, ispace will be modified in the loop. */

iflag:1$
print("Please give floating numerical values to parameters.");
ipar[7]:1$   /* to suppress abundant colnew printout */
break();

/* Example3: Print ltol:[1,3]$ tol:[1d-7,1d-7]$ s:0.2d0$ L:60d0$ followed by exit; */

define(funmake (fsub, cons(var0,var_list)), apply (f, cons(var0,var_list)))$
define(funmake (dfsub, cons(var0,var_list)), apply(df,cons(var0,var_list)))$ 
define(funmake (gsub, cons(i,var_list)), buildq([gval:apply (g, var_list)], gval[i]))$
define(funmake (dgsub, cons(i,var_list)), buildq([dgval:apply(dg,var_list)], row(dgval,i)[1]))$
/* gval and dgval are fully evaluated, [1] converts row to a list */

ipar[7]:1$   /* to suppress colnew printout */

[iflag, fspace, ispace] : colnew_expert(ncomp, m_list, aleft, aright, ZETA,
  ipar, ltol, tol, fixpnt, ispace, fspace, iflag, fsub, dfsub, gsub,
  dgsub, guess)$



/* To get the computed solution at a list xx of points x in [aleft,aright] run
   colnew_appsln. Then zvals is the coresponding list of values of the
   z-vector at these points x. */

/* xx:makelist(0.01*k,k,0,100)$
zvals:colnew_appsln(xx, mstar, fspace, ispace)$
hvals:makelist(zvals[k][1],k,1,101)$
print("Plotting the solution.")$
plot2d([discrete,L*xx,hvals])$  */

Loading /home/michel/.maxima/maxima-init.lisp
Loading /home/michel/.maxima/maxima-init.mac
Maxima 5.48.1 https://maxima.sourceforge.io
using Lisp SBCL 2.5.8
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1) y:1.9d0$

(%i2) aleft : 0.0d0$

(%i3) aright : y$

(%i4) ntol:4$

(%i5) depends([p,f2,f3],t)$

(%i6) lvfs:[t,p,f2,f3]$

(%i7) list_of_eqns: [diff(p,t,2)= -0.5*(1/t+1/(t-1.0)+1/(t-y))*diff(p,t,1)+
((v-f2-f3)/t+f2/(t-1.0)+f3/(t-y))*p , diff(f2,t,1)=0, diff(f3,t,1)=0]$

(%i8) list_of_conds:[(diff(p,t,1)-2.0*p*(v-f2-f3))(0)=0,
(p -1.0)(1)=0, (diff(p,t,1)-2.0*p*f2)(1)=0, (diff(p,t,1)-2.0*p*f3)(y)=0]$

(%i9) list_of_guess:[2.0*t-1.0,1.0,1/(2.0*y-1.0)]$

(%i10) load(colnew)$

(%i11) load("colnew_discov.lisp")$

(%i12) print("Define the interval [aleft,aright] for the differential equations.");
Define the interval [aleft,aright] for the differential equations. 
(%o12) Define the interval [aleft,aright] for the differential equations.
(%i13) if (aleft = 'aleft) then
aleft: float(read("aleft: "));
(%o13)                               false
(%i14) if (aright = 'aright) then
aright: float(read("aright: "));
(%o14)                               false
(%i15) print("Fix the number ntol of checked tolerances");
Fix the number ntol of checked tolerances 
(%o15)             Fix the number ntol of checked tolerances
(%i16) if (ntol = 'ntol) then ntol: read("ntol: ");
(%o16)                               false
(%i17) print("Give the list of variable and functions in equations, like [x,u,v]");
Give the list of variable and functions in equations, like [x,u,v] 
(%o17) Give the list of variable and functions in equations, like [x,u,v]
(%i18) if (lvfs='lvfs) then lvfs :read("lvfs = ");
(%o18)                               false
(%i19) depends(rest(lvfs),first(lvfs))$

(%i20) print("Define the list of the differential equations.");
Define the list of the differential equations. 
(%o20)          Define the list of the differential equations.
(%i21) if (list_of_eqns = 'list_of_eqns) then
list_of_eqns: read("list_of_eqns: ");
(%o21)                               false
(%i22) eqns_to_funs(list_of_eqns)$

(%i23) define(funmake(df,cons(var0,var_list)), jacobian(apply(f,cons(var0,var_list)),var_list))$

(%i24) print("Define the ", mstar, "boundary conditions.");
Define the  4 boundary conditions. 
(%o24)                       boundary conditions.
(%i25) if (list_of_conds = 'list_of_conds) then
list_of_conds: read("list_of_conds: ");
(%o25)                               false
(%i26) conds_to_zconds(list_of_conds, aleft, aright, max_subst_list, m_list)$

(%i27) fixpnt: float(fixpnt)$

(%i28) ZETA:float(ZETA)$

(%i29) define(funmake(dg, var_list), jacobian(apply(g,var_list),var_list))$

(%i30) if (list_of_guess = 'list_of_guess) then
list_of_guess: read("list_of_guess = ");
(%o30)                               false
(%i31) define(guess(var0), build_guess(list_of_guess,m_list,ncomp,var0))$

(%i32) if ntol='ntol then ntol: mstar;
(%o32)                               false
(%i33) build_ipar(m_list, mstar, ncomp, nfixpnt,ntol)$
The array tol is filled by default, you can edit tol 
(%i34) fspace:makelist(0d0,k,1,ipar[5])$

(%i35) ispace:makelist(0,k,1,ipar[6])$

(%i36) print("ncomp = ", ncomp, ", m_list = ", m_list, ", mstar = ", mstar,
", ntol = ", ntol, ", max_subst_list = ", max_subst_list, "");
ncomp =  3 , m_list =  [2, 1, 1] , mstar =  4 , ntol =  4 , max_subst_list =  
                                                             dp
                                          [f3 = z4, f2 = z3, ── = z2, p = z1]  
                                                             dt
(%o36) 
(%i37) print("ZETA: ", ZETA);
ZETA:  [0.0, 1.0, 1.0, 1.9] 
(%o37)                       [0.0, 1.0, 1.0, 1.9]
(%i38) print("ipar: ", ipar);
ipar:  [1, 3, 5, 4, 40000, 2500, 0, 0, 1, 0, 1] 
(%o38)             [1, 3, 5, 4, 40000, 2500, 0, 0, 1, 0, 1]
(%i39) print("ltol: ", ltol);
ltol:  [1, 2, 3, 4] 
(%o39)                           [1, 2, 3, 4]
(%i40) print("tol: ",tol);
tol:  [1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5] 
(%o40)                 [1.0e-5, 1.0e-5, 1.0e-5, 1.0e-5]
(%i41) print("fixpnt: ",fixpnt);
fixpnt:  [1.0] 
(%o41)                               [1.0]
(%i42) print("nfixpnt: ",nfixpnt);
nfixpnt:  1 
(%o42)                                 1
(%i43) iflag:1$

(%i44) print("Please give floating numerical values to parameters.");
Please give floating numerical values to parameters. 
(%o44)       Please give floating numerical values to parameters.
(%i45) ipar[7]:1$

(%i46) break();

Entering a Maxima break point. Type 'exit;' to resume.
_exit;
(%o46)                                []
(%i47) nsteps:100$

(%i48) eigenval: makelist([0.0,0.0,0.0],k,1,nsteps)$

(%i49) v:0.00001d0$

(%i50) for ii from 1 thru nsteps do (
define(funmake (fsub, cons(var0,var_list)), apply (f, cons(var0,var_list))),
define(funmake (dfsub, cons(var0,var_list)), apply(df,cons(var0,var_list))),
define(funmake (gsub, cons(i,var_list)), buildq([gval:apply (g, var_list)], gval[i])),
define(funmake (dgsub, cons(i,var_list)), buildq([dgval:apply(dg,var_list)], row(dgval,i)[1])),
[iflag, fspace, ispace] : colnew_expert(ncomp, m_list, aleft, aright, ZETA,
  ipar, ltol, tol, fixpnt, ispace, fspace, iflag, fsub, dfsub, gsub,
  dgsub, guess),
zval: colnew_appsln([0.1d0], 4, fspace, ispace)[1],
eigenval[ii][1]: zval[3],  /* zval[1] zval[2] are p and p', useless. */
eigenval[ii][2]: zval[4],
eigenval[ii][3]: zval[3]+y*zval[4],
ipar[9]:3,
ipar[3]: ispace[1],
v:v+0.2)$

(%i51) time(%);
(%o51)                            [3.754269]
(%i52) vv: makelist(0.2*k,k,0,nsteps-1)$

(%i53) plot2d([[discrete,vv, makelist(eigenval[k][1],k,1,nsteps)],
  [discrete,vv, makelist(eigenval[k][2],k,1,nsteps)],
  [discrete,vv, makelist(eigenval[k][3],k,1,nsteps)]])$

(%i54)