CoCalc -- MIWsim1DjacobianV1b.sagews

###MIWsim1D###

using lsode with Jacobian

Author: D Cyganski and Bill Page

Created: 2015-05-17

Revised: 2015-06-08

####Implementation of Wiseman Many-Interacting-Worlds simulation for 1D, 1 Particle systems####

single particle with gaussian distribution and nonzero mean velocity, colliding with a potential barrier

The future track of the many worlds position state $Xp(t)=[x_1(t)...x_{N*J*D}(t)]$ is determined in general in MIW by system of $N*J*D$ coupled non-linear second order partial differential equations, where $N$ =number of worlds, $J$ =number of particles, $D$ =number of dimensions. $m_k*\ddot{x}_{kn}(t) = f_k(x_n(t)) + r_{kN}(x_n(t);Xp(t))$ where $k$ indexes the particle and $n$ the world.

Such a system of $N*J*D$ second order non-linear PDEs can always be written as a system of of $2*N*J*D$ first order linear PDEs:

\dot{x}_{kn}(t) = y_{kn}(t) \\ \dot{y}_{kn}(t) = (1/m_k)*(f_k(x_n(t)) + r_{kN}(x_n(t))

For our 1 D case, there will be $2N$ equations.

Setup Sage environment

%sage
# Better version of FriCAS/Axiom mode for SageMathCloud
os.environ['PATH'] = '%s/bin:%s'%(os.environ['HOME'],os.environ['PATH'])
from fricas import fricas
execfile('fricas_md.py')
%default_mode fricas_md

)version
)clear completely
)set output tex on
)set output algebra off
)set message type off

Value = "FriCAS 2014-12-18 compiled at Tue May 5 18:13:01 UTC 2015"

All user variables and function definitions have been cleared.

All )browse facility databases have been cleared.

Internally cached functions and constructors have been cleared.

)clear completely is finished.

Classical Force centered at $x=c$ that repells the oncoming particles

--F(x) == b*(x-c)/(a^3*%pi)*exp(-(x-c)^2/a^2);
F(x) == 9*sech(5*(x-c))^2*tanh(5*(x-c))
F(x)

Compiling function F with type Variable(x) -> Expression(Integer)

$$9 \ {{{\mathrm{sech} \left( {{{5 \ x} -{5 \ c}}} \right)}} ^{2}} \ {\mathrm{tanh} \left( {{{5 \ x} -{5 \ c}}} \right)}$$

complexNormalize htrigs F(x)

$${{{36} \ {{{{e} ^{{{5 \ x} -{5 \ c}}}}} ^{4}}} -{{36} \ {{{{e} ^{{{5 \ x} -{5 \ c}}}}} ^{2}}}} \over {{{{{e} ^{{{5 \ x} -{5 \ c}}}}} ^{6}}+{3 \ {{{{e} ^{{{5 \ x} -{5 \ c}}}}} ^{4}}}+{3 \ {{{{e} ^{{{5 \ x} -{5 \ c}}}}} ^{2}}}+1}$$

%sage
plot(fricas("eval(F(x),[a=0.25,c=5,b=0.7])").sage(),(x,0,10),color="red")

Potential

%sage
plot(fricas("eval(-integrate(F(x),x),[a=0.25,c=5,b=0.7])").sage(),(x,0,10),color="red")

Position of particle in the $i^{th}$ world: $x(i)$

x:=operator('x);

Many Worlds Interaction

Equations 24,25 of HEW paper

R(x,i) == (hbar^2/(4*m))*(sigma(x,i+1)-sigma(x,i))
sigma(x,i) == sif(x,i)^2*(sif(x,i+1)-2*sif(x,i)+sif(x,i-1))
sif(x,i) == 1/(x(i)-x(i-1))

Equations of Motion

xdd(i)==(F(x(i))+R(x,i))/m

The non-constant elements of the Jacobian matrix are formed from derivates with respect to x

jac1(i,j)==realElementary eval(D(eval(xdd(i),x(i+j)=%x),%x),%x=x(i+j))

Sage routine to convert FriCAS expression to Octave function

%sage
# fold lines for octave
def split_lines(fr):
  i=0; cpl=80; f=""
  while i<len(fr):
    if i+cpl>=len(fr):
      f+= "    "+fr[i:i+cpl]+";\n"
      i=i+cpl
    else:
      line=fr[i:i+cpl]
      p=line.rfind('*', 0, cpl)+1
      if p==0:
        p=line.rindex('+', 0, cpl)+1
      f+= "    "+fr[i:i+p]+"...\n"
      i=i+p
  return f[4:len(f)-1]

%sage
# name of function defined by fricas expression ex
def fricas2octave(name,ex):
  fr=fricas(ex).unparsed_input_form()
  # cleanup fricas expresssion
  fr = fr.replace("\n","")
  fr = fr.replace("\r","")
  fr = fr.replace(" -","-")
  fr = fr.replace("+-","-")
  # boundry conditions
  fr = fr.replace("x(i-2)","xl2")
  fr = fr.replace("x(i-1)","xl1")
  fr = fr.replace("x(i+1)","xr1")
  fr = fr.replace("x(i+2)","xr2")
  f = """\
function result=%s(x,i)
  global l2 l1 p1 p2;
  global N m hbar a b c;
  if (i-2<1)
    xl2=l2;
  else
    xl2=x(i-2);
  endif;
  if (i-1<1)
    xl1=l1;
  else
    xl1=x(i-1);
  endif;
  if (i+2>N)
    xr2=p2;
  else
    xr2=x(i+2);
  endif;
  if (i+1>N)
    xr1=p1;
  else
    xr1=x(i+1);
  endif;
  result=%s
endfunction
"""%(name,split_lines(fr))
  ff = open('%s.m'%name, 'w')
  ff.write(f)
  ff.close()
  return f

%sage
# Use a better version of the octave interface
from octave import octave
%default_mode octave

Generate initial distribution of particle states, including position,x (X0), and velocity, xdot (Y0)

clear all;
version

ans = 3.8.2

For a simple test, distribute 101 values equally over a window and set all velocities to zero

#xvalues=[-1:.05:1]'; #Prob2
#xvalues=[-1:.025:1]'; #Prob3
#xvalues=[-0.5:.025:0.5]'; #Prob4

Use inverse error function to generate a Gaussian distribution of trajectories

#xuniform=[-0.95:0.1:0.95]'; %For the simple test
#xuniform=[-0.95:0.022|:0.95]';
xuniform=[-0.99:0.15:0.99]';
#xvalues=erfcinv(xuniform);
xvalues=erfinv(xuniform);
global N; %Number of worlds
N = size(xvalues,1)
X0 = xvalues;
#Y0 = zeros(N,1); %zero initial velocities
#Y0 = normrnd(1, 0.2, N, 1);
Y0 = ones(N,1); %non-zero initial velocities

N = 14

Plot the input trajectory density

Pin=(N*(xvalues(2:N,1)-xvalues(1:(N-1),1))).^(-1);

%sage list_plot(matrix(RR,octave("[xvalues(2:N),Pin]")),plotjoined=true)

Constants

global m hbar;
m = 1; %mass of particle
hbar = 0.1; %Planck constant %.02 works but nonguassian after 500 iterations
global a b c;
a = 0.25;
b = 0.7;
c = 5;

Boundry conditions

global l2 l1 p1 p2;
l2 = -100000;
l1 =  -10000;
p1 =   10000;
p2 =  100000;

%sage
octave.eval("global F = @(x) "+fricas("F('x)").unparsed_input_form().replace("pi()","pi"))

''

Plot the force

Xrange = 0:0.05:10;

%sage list_plot(matrix(RR,octave("[Xrange',arrayfun(F,Xrange')]")),plotjoined=true)

%sage print fricas2octave('R2','R(x,i)')

function result=R2(x,i)
  global l2 l1 p1 p2;
  global N m hbar a b c;
  if (i-2<1)
    xl2=l2;
  else
    xl2=x(i-2);
  endif;
  if (i-1<1)
    xl1=l1;
  else
    xl1=x(i-1);
  endif;
  if (i+2>N)
    xr2=p2;
  else
    xr2=x(i+2);
  endif;
  if (i+1>N)
    xr1=p1;
  else
    xr1=x(i+1);
  endif;
  result=(((-1*hbar^2*x(i)+(3*hbar^2*xl1-2*hbar^2*xl2))*xr1^3+(3*hbar^2*x(i)^2+(-10*...
    hbar^2*xl1+7*hbar^2*xl2)*x(i)+(hbar^2*xl1^2-1*hbar^2*xl2*xl1))*xr1^2+(-3*hbar^2*...
    x(i)^3+(12*hbar^2*xl1-9*hbar^2*xl2)*x(i)^2+(-4*hbar^2*xl1^2+4*hbar^2*xl2*xl1)*...
    x(i)+(hbar^2*xl1^3-1*hbar^2*xl2*xl1^2))*xr1+(hbar^2*x(i)^4+(-7*hbar^2*xl1+6*...
    hbar^2*xl2)*x(i)^3+(9*hbar^2*xl1^2-9*hbar^2*xl2*xl1)*x(i)^2+(-7*hbar^2*xl1^3+7*...
    hbar^2*xl2*xl1^2)*x(i)+(2*hbar^2*xl1^4-2*hbar^2*xl2*xl1^3)))*xr2+((hbar^2*...
    x(i)+(-3*hbar^2*xl1+2*hbar^2*xl2))*xr1^4+(-3*hbar^2*x(i)^2+(10*hbar^2*xl1-7*...
    hbar^2*xl2)*x(i)+(-1*hbar^2*xl1^2+hbar^2*xl2*xl1))*xr1^3+(3*hbar^2*x(i)^3+(-12*...
    hbar^2*xl1+9*hbar^2*xl2)*x(i)^2+(4*hbar^2*xl1^2-4*hbar^2*xl2*xl1)*x(i)+(-1*...
    hbar^2*xl1^3+hbar^2*xl2*xl1^2))*xr1^2+(-1*hbar^2*x(i)^4+(8*hbar^2*xl1-7*hbar^2*...
    xl2)*x(i)^3+(-12*hbar^2*xl1^2+12*hbar^2*xl2*xl1)*x(i)^2+(10*hbar^2*xl1^3-10*...
    hbar^2*xl2*xl1^2)*x(i)+(-3*hbar^2*xl1^4+3*hbar^2*xl2*xl1^3))*xr1+((-1*hbar^2*...
    xl1+hbar^2*xl2)*x(i)^4+(3*hbar^2*xl1^2-3*hbar^2*xl2*xl1)*x(i)^3+(-3*hbar^2*...
    xl1^3+3*hbar^2*xl2*xl1^2)*x(i)^2+(hbar^2*xl1^4-1*hbar^2*xl2*xl1^3)*...
    x(i))))/((((4*m*xl1-4*m*xl2)*x(i)^3+(-12*m*xl1^2+12*m*xl2*xl1)*x(i)^2+(12*m*...
    xl1^3-12*m*xl2*xl1^2)*x(i)+(-4*m*xl1^4+4*m*xl2*xl1^3))*xr1^3+((-12*m*xl1+12*m*...
    xl2)*x(i)^4+(36*m*xl1^2-36*m*xl2*xl1)*x(i)^3+(-36*m*xl1^3+36*m*xl2*xl1^2)*...
    x(i)^2+(12*m*xl1^4-12*m*xl2*xl1^3)*x(i))*xr1^2+((12*m*xl1-12*m*xl2)*x(i)^5+(-36*...
    m*xl1^2+36*m*xl2*xl1)*x(i)^4+(36*m*xl1^3-36*m*xl2*xl1^2)*x(i)^3+(-12*m*xl1^4+12*...
    m*xl2*xl1^3)*x(i)^2)*xr1+((-4*m*xl1+4*m*xl2)*x(i)^6+(12*m*xl1^2-12*m*xl2*xl1)*...
    x(i)^5+(-12*m*xl1^3+12*m*xl2*xl1^2)*x(i)^4+(4*m*xl1^4-4*m*xl2*xl1^3)*x(i)^3))*...
    xr2+(((-4*m*xl1+4*m*xl2)*x(i)^3+(12*m*xl1^2-12*m*xl2*xl1)*x(i)^2+(-12*m*...
    xl1^3+12*m*xl2*xl1^2)*x(i)+(4*m*xl1^4-4*m*xl2*xl1^3))*xr1^4+((12*m*xl1-12*m*...
    xl2)*x(i)^4+(-36*m*xl1^2+36*m*xl2*xl1)*x(i)^3+(36*m*xl1^3-36*m*xl2*xl1^2)*...
    x(i)^2+(-12*m*xl1^4+12*m*xl2*xl1^3)*x(i))*xr1^3+((-12*m*xl1+12*m*xl2)*...
    x(i)^5+(36*m*xl1^2-36*m*xl2*xl1)*x(i)^4+(-36*m*xl1^3+36*m*xl2*xl1^2)*x(i)^3+(12*...
    m*xl1^4-12*m*xl2*xl1^3)*x(i)^2)*xr1^2+((4*m*xl1-4*m*xl2)*x(i)^6+(-12*m*xl1^2+12*...
    m*xl2*xl1)*x(i)^5+(12*m*xl1^3-12*m*xl2*xl1^2)*x(i)^4+(-4*m*xl1^4+4*m*xl2*xl1^3)*...
    x(i)^3)*xr1));
endfunction

R2(X0,3)

ans = 0.00280219

function r = R(x,i)
  global m hbar;
  r = (hbar^2/(4*m))*(sigma(x,i+1)-sigma(x,i));
endfunction;

function s=sigma(x,i)
 s = sif(x,i)^2*(sif(x,i+1)-2*sif(x,i)+sif(x,i-1));
endfunction;

function den = sif(x,i)
  global N
  if (i>N)||(i<2)
    den=0;
  else
    den=1/(x(i)-x(i-1));
  endif
endfunction;

R2(X0,3)

ans = 0.00280219

%sage
X=map(eval,str(octave('X0')).splitlines()[2:-1]);X

[-1.82139, -0.993536, -0.71787, -0.522444, -0.360676, -0.216008, -0.0799303, 0.0532238, 0.18831, 0.330713, 0.488122, 0.674697, 0.926719, 1.45222]

%sage
print fricas.eval("X:=%s;"%(fricas(X).name()))

%fricas_md
eval(R(X,3),[hbar=0.1,m=1])

Compiling function sif with type (List(Float),PositiveInteger) -> Float

Compiling function sif with type (List(Float),Integer) -> Float

Compiling function sigma with type (List(Float),PositiveInteger) -> Float

Compiling function R with type (List(Float),PositiveInteger) -> Fraction(Polynomial(Float))

0.0027954071\_1864352625\_97

Equations of motion in numeric form

function xdot = deq(x,t)
    global F N;
    for i=1:N
      xdot(i) = x(N+i);
      xdot(N+i) = F(x(i))+R(x,i);
    endfor
endfunction

Functions to compute the Jacobian matrix

%sage print fricas2octave('Jm2','jac1(i,-2)')

function result=Jm2(x,i)
  global l2 l1 p1 p2;
  global N m hbar a b c;
  if (i-2<1)
    xl2=l2;
  else
    xl2=x(i-2);
  endif;
  if (i-1<1)
    xl1=l1;
  else
    xl1=x(i-1);
  endif;
  if (i+2>N)
    xr2=p2;
  else
    xr2=x(i+2);
  endif;
  if (i+1>N)
    xr1=p1;
  else
    xr1=x(i+1);
  endif;
  result=(-1*hbar^2)/((4*m^2*xl1^2-8*m^2*xl2*xl1+4*m^2*xl2^2)*x(i)^2+(-8*m^2*xl1^3+16*...
    m^2*xl2*xl1^2-8*m^2*xl2^2*xl1)*x(i)+(4*m^2*xl1^4-8*m^2*xl2*xl1^3+4*m^2*xl2^2*...
    xl1^2));
endfunction

%sage print fricas2octave('Jm1','jac1(i,-1)')

function result=Jm1(x,i)
  global l2 l1 p1 p2;
  global N m hbar a b c;
  if (i-2<1)
    xl2=l2;
  else
    xl2=x(i-2);
  endif;
  if (i-1<1)
    xl1=l1;
  else
    xl1=x(i-1);
  endif;
  if (i+2>N)
    xr2=p2;
  else
    xr2=x(i+2);
  endif;
  if (i+1>N)
    xr1=p1;
  else
    xr1=x(i+1);
  endif;
  result=((hbar^2*x(i)^2+(-4*hbar^2*xl1+2*hbar^2*xl2)*x(i)+(9*hbar^2*xl1^2-14*hbar^2*xl2*...
    xl1+6*hbar^2*xl2^2))*xr1^2+(-2*hbar^2*x(i)^3+(8*hbar^2*xl1-4*hbar^2*xl2)*...
    x(i)^2+(-20*hbar^2*xl1^2+32*hbar^2*xl2*xl1-14*hbar^2*xl2^2)*x(i)+(2*hbar^2*...
    xl1^3-4*hbar^2*xl2*xl1^2+2*hbar^2*xl2^2*xl1))*xr1+(hbar^2*x(i)^4+(-4*hbar^2*...
    xl1+2*hbar^2*xl2)*x(i)^3+(12*hbar^2*xl1^2-20*hbar^2*xl2*xl1+9*hbar^2*xl2^2)*...
    x(i)^2+(-4*hbar^2*xl1^3+8*hbar^2*xl2*xl1^2-4*hbar^2*xl2^2*xl1)*x(i)+(hbar^2*...
    xl1^4-2*hbar^2*xl2*xl1^3+hbar^2*xl2^2*xl1^2)))/(((4*m^2*xl1^2-8*m^2*xl2*xl1+4*...
    m^2*xl2^2)*x(i)^4+(-16*m^2*xl1^3+32*m^2*xl2*xl1^2-16*m^2*xl2^2*xl1)*x(i)^3+(24*...
    m^2*xl1^4-48*m^2*xl2*xl1^3+24*m^2*xl2^2*xl1^2)*x(i)^2+(-16*m^2*xl1^5+32*m^2*xl2*...
    xl1^4-16*m^2*xl2^2*xl1^3)*x(i)+(4*m^2*xl1^6-8*m^2*xl2*xl1^5+4*m^2*xl2^2*xl1^4))*...
    xr1^2+((-8*m^2*xl1^2+16*m^2*xl2*xl1-8*m^2*xl2^2)*x(i)^5+(32*m^2*xl1^3-64*m^2*...
    xl2*xl1^2+32*m^2*xl2^2*xl1)*x(i)^4+(-48*m^2*xl1^4+96*m^2*xl2*xl1^3-48*m^2*xl2^2*...
    xl1^2)*x(i)^3+(32*m^2*xl1^5-64*m^2*xl2*xl1^4+32*m^2*xl2^2*xl1^3)*x(i)^2+(-8*m^2*...
    xl1^6+16*m^2*xl2*xl1^5-8*m^2*xl2^2*xl1^4)*x(i))*xr1+((4*m^2*xl1^2-8*m^2*xl2*...
    xl1+4*m^2*xl2^2)*x(i)^6+(-16*m^2*xl1^3+32*m^2*xl2*xl1^2-16*m^2*xl2^2*xl1)*...
    x(i)^5+(24*m^2*xl1^4-48*m^2*xl2*xl1^3+24*m^2*xl2^2*xl1^2)*x(i)^4+(-16*m^2*...
    xl1^5+32*m^2*xl2*xl1^4-16*m^2*xl2^2*xl1^3)*x(i)^3+(4*m^2*xl1^6-8*m^2*xl2*...
    xl1^5+4*m^2*xl2^2*xl1^4)*x(i)^2));
endfunction

%sage print fricas2octave('J0','jac1(i,0)')

function result=J0(x,i)
  global l2 l1 p1 p2;
  global N m hbar a b c;
  if (i-2<1)
    xl2=l2;
  else
    xl2=x(i-2);
  endif;
  if (i-1<1)
    xl1=l1;
  else
    xl1=x(i-1);
  endif;
  if (i+2>N)
    xr2=p2;
  else
    xr2=x(i+2);
  endif;
  if (i+1>N)
    xr1=p1;
  else
    xr1=x(i+1);
  endif;
  result=((((hbar^2*x(i)+(-4*hbar^2*xl1+3*hbar^2*xl2))*xr1^4+(-4*hbar^2*x(i)^2+(17*...
    hbar^2*xl1-13*hbar^2*xl2)*x(i)+(-1*hbar^2*xl1^2+hbar^2*xl2*xl1))*xr1^3+(6*...
    hbar^2*x(i)^3+(-28*hbar^2*xl1+22*hbar^2*xl2)*x(i)^2+(5*hbar^2*xl1^2-5*hbar^2*...
    xl2*xl1)*x(i)+(-1*hbar^2*xl1^3+hbar^2*xl2*xl1^2))*xr1^2+(-4*hbar^2*x(i)^4+(22*...
    hbar^2*xl1-18*hbar^2*xl2)*x(i)^3+(-10*hbar^2*xl1^2+10*hbar^2*xl2*xl1)*x(i)^2+(5*...
    hbar^2*xl1^3-5*hbar^2*xl2*xl1^2)*x(i)+(-1*hbar^2*xl1^4+hbar^2*xl2*xl1^3))*...
    xr1+(hbar^2*x(i)^5+(-10*hbar^2*xl1+9*hbar^2*xl2)*x(i)^4+(18*hbar^2*xl1^2-18*...
    hbar^2*xl2*xl1)*x(i)^3+(-22*hbar^2*xl1^3+22*hbar^2*xl2*xl1^2)*x(i)^2+(13*hbar^2*...
    xl1^4-13*hbar^2*xl2*xl1^3)*x(i)+(-3*hbar^2*xl1^5+3*hbar^2*xl2*xl1^4)))*xr2+((-1*...
    hbar^2*x(i)+(4*hbar^2*xl1-3*hbar^2*xl2))*xr1^5+(4*hbar^2*x(i)^2+(-17*hbar^2*...
    xl1+13*hbar^2*xl2)*x(i)+(hbar^2*xl1^2-1*hbar^2*xl2*xl1))*xr1^4+(-6*hbar^2*...
    x(i)^3+(28*hbar^2*xl1-22*hbar^2*xl2)*x(i)^2+(-5*hbar^2*xl1^2+5*hbar^2*xl2*xl1)*...
    x(i)+(hbar^2*xl1^3-1*hbar^2*xl2*xl1^2))*xr1^3+(4*hbar^2*x(i)^4+(-22*hbar^2*...
    xl1+18*hbar^2*xl2)*x(i)^3+(10*hbar^2*xl1^2-10*hbar^2*xl2*xl1)*x(i)^2+(-5*hbar^2*...
    xl1^3+5*hbar^2*xl2*xl1^2)*x(i)+(hbar^2*xl1^4-1*hbar^2*xl2*xl1^3))*xr1^2+(-1*...
    hbar^2*x(i)^5+(11*hbar^2*xl1-10*hbar^2*xl2)*x(i)^4+(-22*hbar^2*xl1^2+22*hbar^2*...
    xl2*xl1)*x(i)^3+(28*hbar^2*xl1^3-28*hbar^2*xl2*xl1^2)*x(i)^2+(-17*hbar^2*...
    xl1^4+17*hbar^2*xl2*xl1^3)*x(i)+(4*hbar^2*xl1^5-4*hbar^2*xl2*xl1^4))*xr1+((-1*...
    hbar^2*xl1+hbar^2*xl2)*x(i)^5+(4*hbar^2*xl1^2-4*hbar^2*xl2*xl1)*x(i)^4+(-6*...
    hbar^2*xl1^3+6*hbar^2*xl2*xl1^2)*x(i)^3+(4*hbar^2*xl1^4-4*hbar^2*xl2*xl1^3)*...
    x(i)^2+(-1*hbar^2*xl1^5+hbar^2*xl2*xl1^4)*x(i))))*exp(5*x(i)-5*c)^8+((((-720*m*...
    xl1+720*m*xl2)*x(i)^4+(2880*m*xl1^2-2880*m*xl2*xl1)*x(i)^3+(-4320*m*xl1^3+4320*...
    m*xl2*xl1^2)*x(i)^2+(2880*m*xl1^4-2880*m*xl2*xl1^3+4*hbar^2)*x(i)+(-720*m*...
    xl1^5+720*m*xl2*xl1^4-16*hbar^2*xl1+12*hbar^2*xl2))*xr1^4+((2880*m*xl1-2880*m*...
    xl2)*x(i)^5+(-11520*m*xl1^2+11520*m*xl2*xl1)*x(i)^4+(17280*m*xl1^3-17280*m*xl2*...
    xl1^2)*x(i)^3+(-11520*m*xl1^4+11520*m*xl2*xl1^3-16*hbar^2)*x(i)^2+(2880*m*...
    xl1^5-2880*m*xl2*xl1^4+68*hbar^2*xl1-52*hbar^2*xl2)*x(i)+(-4*hbar^2*xl1^2+4*...
    hbar^2*xl2*xl1))*xr1^3+((-4320*m*xl1+4320*m*xl2)*x(i)^6+(17280*m*xl1^2-17280*m*...
    xl2*xl1)*x(i)^5+(-25920*m*xl1^3+25920*m*xl2*xl1^2)*x(i)^4+(17280*m*xl1^4-17280*...
    m*xl2*xl1^3+24*hbar^2)*x(i)^3+(-4320*m*xl1^5+4320*m*xl2*xl1^4-112*hbar^2*xl1+88*...
    hbar^2*xl2)*x(i)^2+(20*hbar^2*xl1^2-20*hbar^2*xl2*xl1)*x(i)+(-4*hbar^2*xl1^3+4*...
    hbar^2*xl2*xl1^2))*xr1^2+((2880*m*xl1-2880*m*xl2)*x(i)^7+(-11520*m*xl1^2+11520*...
    m*xl2*xl1)*x(i)^6+(17280*m*xl1^3-17280*m*xl2*xl1^2)*x(i)^5+(-11520*m*...
    xl1^4+11520*m*xl2*xl1^3-16*hbar^2)*x(i)^4+(2880*m*xl1^5-2880*m*xl2*xl1^4+88*...
    hbar^2*xl1-72*hbar^2*xl2)*x(i)^3+(-40*hbar^2*xl1^2+40*hbar^2*xl2*xl1)*...
    x(i)^2+(20*hbar^2*xl1^3-20*hbar^2*xl2*xl1^2)*x(i)+(-4*hbar^2*xl1^4+4*hbar^2*xl2*...
    xl1^3))*xr1+((-720*m*xl1+720*m*xl2)*x(i)^8+(2880*m*xl1^2-2880*m*xl2*xl1)*...
    x(i)^7+(-4320*m*xl1^3+4320*m*xl2*xl1^2)*x(i)^6+(2880*m*xl1^4-2880*m*xl2*xl1^3+4*...
    hbar^2)*x(i)^5+(-720*m*xl1^5+720*m*xl2*xl1^4-40*hbar^2*xl1+36*hbar^2*xl2)*...
    x(i)^4+(72*hbar^2*xl1^2-72*hbar^2*xl2*xl1)*x(i)^3+(-88*hbar^2*xl1^3+88*hbar^2*...
    xl2*xl1^2)*x(i)^2+(52*hbar^2*xl1^4-52*hbar^2*xl2*xl1^3)*x(i)+(-12*hbar^2*...
    xl1^5+12*hbar^2*xl2*xl1^4)))*xr2+(((720*m*xl1-720*m*xl2)*x(i)^4+(-2880*m*...
    xl1^2+2880*m*xl2*xl1)*x(i)^3+(4320*m*xl1^3-4320*m*xl2*xl1^2)*x(i)^2+(-2880*m*...
    xl1^4+2880*m*xl2*xl1^3-4*hbar^2)*x(i)+(720*m*xl1^5-720*m*xl2*xl1^4+16*hbar^2*...
    xl1-12*hbar^2*xl2))*xr1^5+((-2880*m*xl1+2880*m*xl2)*x(i)^5+(11520*m*xl1^2-11520*...
    m*xl2*xl1)*x(i)^4+(-17280*m*xl1^3+17280*m*xl2*xl1^2)*x(i)^3+(11520*m*...
    xl1^4-11520*m*xl2*xl1^3+16*hbar^2)*x(i)^2+(-2880*m*xl1^5+2880*m*xl2*xl1^4-68*...
    hbar^2*xl1+52*hbar^2*xl2)*x(i)+(4*hbar^2*xl1^2-4*hbar^2*xl2*xl1))*xr1^4+((4320*...
    m*xl1-4320*m*xl2)*x(i)^6+(-17280*m*xl1^2+17280*m*xl2*xl1)*x(i)^5+(25920*m*...
    xl1^3-25920*m*xl2*xl1^2)*x(i)^4+(-17280*m*xl1^4+17280*m*xl2*xl1^3-24*hbar^2)*...
    x(i)^3+(4320*m*xl1^5-4320*m*xl2*xl1^4+112*hbar^2*xl1-88*hbar^2*xl2)*x(i)^2+(-20*...
    hbar^2*xl1^2+20*hbar^2*xl2*xl1)*x(i)+(4*hbar^2*xl1^3-4*hbar^2*xl2*xl1^2))*...
    xr1^3+((-2880*m*xl1+2880*m*xl2)*x(i)^7+(11520*m*xl1^2-11520*m*xl2*xl1)*...
    x(i)^6+(-17280*m*xl1^3+17280*m*xl2*xl1^2)*x(i)^5+(11520*m*xl1^4-11520*m*xl2*...
    xl1^3+16*hbar^2)*x(i)^4+(-2880*m*xl1^5+2880*m*xl2*xl1^4-88*hbar^2*xl1+72*hbar^2*...
    xl2)*x(i)^3+(40*hbar^2*xl1^2-40*hbar^2*xl2*xl1)*x(i)^2+(-20*hbar^2*xl1^3+20*...
    hbar^2*xl2*xl1^2)*x(i)+(4*hbar^2*xl1^4-4*hbar^2*xl2*xl1^3))*xr1^2+((720*m*...
    xl1-720*m*xl2)*x(i)^8+(-2880*m*xl1^2+2880*m*xl2*xl1)*x(i)^7+(4320*m*xl1^3-4320*...
    m*xl2*xl1^2)*x(i)^6+(-2880*m*xl1^4+2880*m*xl2*xl1^3-4*hbar^2)*x(i)^5+(720*m*...
    xl1^5-720*m*xl2*xl1^4+44*hbar^2*xl1-40*hbar^2*xl2)*x(i)^4+(-88*hbar^2*xl1^2+88*...
    hbar^2*xl2*xl1)*x(i)^3+(112*hbar^2*xl1^3-112*hbar^2*xl2*xl1^2)*x(i)^2+(-68*...
    hbar^2*xl1^4+68*hbar^2*xl2*xl1^3)*x(i)+(16*hbar^2*xl1^5-16*hbar^2*xl2*xl1^4))*...
    xr1+((-4*hbar^2*xl1+4*hbar^2*xl2)*x(i)^5+(16*hbar^2*xl1^2-16*hbar^2*xl2*xl1)*...
    x(i)^4+(-24*hbar^2*xl1^3+24*hbar^2*xl2*xl1^2)*x(i)^3+(16*hbar^2*xl1^4-16*hbar^2*...
    xl2*xl1^3)*x(i)^2+(-4*hbar^2*xl1^5+4*hbar^2*xl2*xl1^4)*x(i))))*exp(5*x(i)-5*...
    c)^6+((((2880*m*xl1-2880*m*xl2)*x(i)^4+(-11520*m*xl1^2+11520*m*xl2*xl1)*...
    x(i)^3+(17280*m*xl1^3-17280*m*xl2*xl1^2)*x(i)^2+(-11520*m*xl1^4+11520*m*xl2*...
    xl1^3+6*hbar^2)*x(i)+(2880*m*xl1^5-2880*m*xl2*xl1^4-24*hbar^2*xl1+18*hbar^2*...
    xl2))*xr1^4+((-11520*m*xl1+11520*m*xl2)*x(i)^5+(46080*m*xl1^2-46080*m*xl2*xl1)*...
    x(i)^4+(-69120*m*xl1^3+69120*m*xl2*xl1^2)*x(i)^3+(46080*m*xl1^4-46080*m*xl2*...
    xl1^3-24*hbar^2)*x(i)^2+(-11520*m*xl1^5+11520*m*xl2*xl1^4+102*hbar^2*xl1-78*...
    hbar^2*xl2)*x(i)+(-6*hbar^2*xl1^2+6*hbar^2*xl2*xl1))*xr1^3+((17280*m*xl1-17280*...
    m*xl2)*x(i)^6+(-69120*m*xl1^2+69120*m*xl2*xl1)*x(i)^5+(103680*m*xl1^3-103680*m*...
    xl2*xl1^2)*x(i)^4+(-69120*m*xl1^4+69120*m*xl2*xl1^3+36*hbar^2)*x(i)^3+(17280*m*...
    xl1^5-17280*m*xl2*xl1^4-168*hbar^2*xl1+132*hbar^2*xl2)*x(i)^2+(30*hbar^2*...
    xl1^2-30*hbar^2*xl2*xl1)*x(i)+(-6*hbar^2*xl1^3+6*hbar^2*xl2*xl1^2))*...
    xr1^2+((-11520*m*xl1+11520*m*xl2)*x(i)^7+(46080*m*xl1^2-46080*m*xl2*xl1)*...
    x(i)^6+(-69120*m*xl1^3+69120*m*xl2*xl1^2)*x(i)^5+(46080*m*xl1^4-46080*m*xl2*...
    xl1^3-24*hbar^2)*x(i)^4+(-11520*m*xl1^5+11520*m*xl2*xl1^4+132*hbar^2*xl1-108*...
    hbar^2*xl2)*x(i)^3+(-60*hbar^2*xl1^2+60*hbar^2*xl2*xl1)*x(i)^2+(30*hbar^2*...
    xl1^3-30*hbar^2*xl2*xl1^2)*x(i)+(-6*hbar^2*xl1^4+6*hbar^2*xl2*xl1^3))*...
    xr1+((2880*m*xl1-2880*m*xl2)*x(i)^8+(-11520*m*xl1^2+11520*m*xl2*xl1)*...
    x(i)^7+(17280*m*xl1^3-17280*m*xl2*xl1^2)*x(i)^6+(-11520*m*xl1^4+11520*m*xl2*...
    xl1^3+6*hbar^2)*x(i)^5+(2880*m*xl1^5-2880*m*xl2*xl1^4-60*hbar^2*xl1+54*hbar^2*...
    xl2)*x(i)^4+(108*hbar^2*xl1^2-108*hbar^2*xl2*xl1)*x(i)^3+(-132*hbar^2*xl1^3+132*...
    hbar^2*xl2*xl1^2)*x(i)^2+(78*hbar^2*xl1^4-78*hbar^2*xl2*xl1^3)*x(i)+(-18*hbar^2*...
    xl1^5+18*hbar^2*xl2*xl1^4)))*xr2+(((-2880*m*xl1+2880*m*xl2)*x(i)^4+(11520*m*...
    xl1^2-11520*m*xl2*xl1)*x(i)^3+(-17280*m*xl1^3+17280*m*xl2*xl1^2)*x(i)^2+(11520*...
    m*xl1^4-11520*m*xl2*xl1^3-6*hbar^2)*x(i)+(-2880*m*xl1^5+2880*m*xl2*xl1^4+24*...
    hbar^2*xl1-18*hbar^2*xl2))*xr1^5+((11520*m*xl1-11520*m*xl2)*x(i)^5+(-46080*m*...
    xl1^2+46080*m*xl2*xl1)*x(i)^4+(69120*m*xl1^3-69120*m*xl2*xl1^2)*x(i)^3+(-46080*...
    m*xl1^4+46080*m*xl2*xl1^3+24*hbar^2)*x(i)^2+(11520*m*xl1^5-11520*m*xl2*...
    xl1^4-102*hbar^2*xl1+78*hbar^2*xl2)*x(i)+(6*hbar^2*xl1^2-6*hbar^2*xl2*xl1))*...
    xr1^4+((-17280*m*xl1+17280*m*xl2)*x(i)^6+(69120*m*xl1^2-69120*m*xl2*xl1)*...
    x(i)^5+(-103680*m*xl1^3+103680*m*xl2*xl1^2)*x(i)^4+(69120*m*xl1^4-69120*m*xl2*...
    xl1^3-36*hbar^2)*x(i)^3+(-17280*m*xl1^5+17280*m*xl2*xl1^4+168*hbar^2*xl1-132*...
    hbar^2*xl2)*x(i)^2+(-30*hbar^2*xl1^2+30*hbar^2*xl2*xl1)*x(i)+(6*hbar^2*xl1^3-6*...
    hbar^2*xl2*xl1^2))*xr1^3+((11520*m*xl1-11520*m*xl2)*x(i)^7+(-46080*m*...
    xl1^2+46080*m*xl2*xl1)*x(i)^6+(69120*m*xl1^3-69120*m*xl2*xl1^2)*x(i)^5+(-46080*...
    m*xl1^4+46080*m*xl2*xl1^3+24*hbar^2)*x(i)^4+(11520*m*xl1^5-11520*m*xl2*...
    xl1^4-132*hbar^2*xl1+108*hbar^2*xl2)*x(i)^3+(60*hbar^2*xl1^2-60*hbar^2*xl2*xl1)*...
    x(i)^2+(-30*hbar^2*xl1^3+30*hbar^2*xl2*xl1^2)*x(i)+(6*hbar^2*xl1^4-6*hbar^2*xl2*...
    xl1^3))*xr1^2+((-2880*m*xl1+2880*m*xl2)*x(i)^8+(11520*m*xl1^2-11520*m*xl2*xl1)*...
    x(i)^7+(-17280*m*xl1^3+17280*m*xl2*xl1^2)*x(i)^6+(11520*m*xl1^4-11520*m*xl2*...
    xl1^3-6*hbar^2)*x(i)^5+(-2880*m*xl1^5+2880*m*xl2*xl1^4+66*hbar^2*xl1-60*hbar^2*...
    xl2)*x(i)^4+(-132*hbar^2*xl1^2+132*hbar^2*xl2*xl1)*x(i)^3+(168*hbar^2*xl1^3-168*...
    hbar^2*xl2*xl1^2)*x(i)^2+(-102*hbar^2*xl1^4+102*hbar^2*xl2*xl1^3)*x(i)+(24*...
    hbar^2*xl1^5-24*hbar^2*xl2*xl1^4))*xr1+((-6*hbar^2*xl1+6*hbar^2*xl2)*x(i)^5+(24*...
    hbar^2*xl1^2-24*hbar^2*xl2*xl1)*x(i)^4+(-36*hbar^2*xl1^3+36*hbar^2*xl2*xl1^2)*...
    x(i)^3+(24*hbar^2*xl1^4-24*hbar^2*xl2*xl1^3)*x(i)^2+(-6*hbar^2*xl1^5+6*hbar^2*...
    xl2*xl1^4)*x(i))))*exp(5*x(i)-5*c)^4+((((-720*m*xl1+720*m*xl2)*x(i)^4+(2880*m*...
    xl1^2-2880*m*xl2*xl1)*x(i)^3+(-4320*m*xl1^3+4320*m*xl2*xl1^2)*x(i)^2+(2880*m*...
    xl1^4-2880*m*xl2*xl1^3+4*hbar^2)*x(i)+(-720*m*xl1^5+720*m*xl2*xl1^4-16*hbar^2*...
    xl1+12*hbar^2*xl2))*xr1^4+((2880*m*xl1-2880*m*xl2)*x(i)^5+(-11520*m*xl1^2+11520*...
    m*xl2*xl1)*x(i)^4+(17280*m*xl1^3-17280*m*xl2*xl1^2)*x(i)^3+(-11520*m*...
    xl1^4+11520*m*xl2*xl1^3-16*hbar^2)*x(i)^2+(2880*m*xl1^5-2880*m*xl2*xl1^4+68*...
    hbar^2*xl1-52*hbar^2*xl2)*x(i)+(-4*hbar^2*xl1^2+4*hbar^2*xl2*xl1))*...
    xr1^3+((-4320*m*xl1+4320*m*xl2)*x(i)^6+(17280*m*xl1^2-17280*m*xl2*xl1)*...
    x(i)^5+(-25920*m*xl1^3+25920*m*xl2*xl1^2)*x(i)^4+(17280*m*xl1^4-17280*m*xl2*...
    xl1^3+24*hbar^2)*x(i)^3+(-4320*m*xl1^5+4320*m*xl2*xl1^4-112*hbar^2*xl1+88*...
    hbar^2*xl2)*x(i)^2+(20*hbar^2*xl1^2-20*hbar^2*xl2*xl1)*x(i)+(-4*hbar^2*xl1^3+4*...
    hbar^2*xl2*xl1^2))*xr1^2+((2880*m*xl1-2880*m*xl2)*x(i)^7+(-11520*m*xl1^2+11520*...
    m*xl2*xl1)*x(i)^6+(17280*m*xl1^3-17280*m*xl2*xl1^2)*x(i)^5+(-11520*m*...
    xl1^4+11520*m*xl2*xl1^3-16*hbar^2)*x(i)^4+(2880*m*xl1^5-2880*m*xl2*xl1^4+88*...
    hbar^2*xl1-72*hbar^2*xl2)*x(i)^3+(-40*hbar^2*xl1^2+40*hbar^2*xl2*xl1)*...
    x(i)^2+(20*hbar^2*xl1^3-20*hbar^2*xl2*xl1^2)*x(i)+(-4*hbar^2*xl1^4+4*hbar^2*xl2*...
    xl1^3))*xr1+((-720*m*xl1+720*m*xl2)*x(i)^8+(2880*m*xl1^2-2880*m*xl2*xl1)*...
    x(i)^7+(-4320*m*xl1^3+4320*m*xl2*xl1^2)*x(i)^6+(2880*m*xl1^4-2880*m*xl2*xl1^3+4*...
    hbar^2)*x(i)^5+(-720*m*xl1^5+720*m*xl2*xl1^4-40*hbar^2*xl1+36*hbar^2*xl2)*...
    x(i)^4+(72*hbar^2*xl1^2-72*hbar^2*xl2*xl1)*x(i)^3+(-88*hbar^2*xl1^3+88*hbar^2*...
    xl2*xl1^2)*x(i)^2+(52*hbar^2*xl1^4-52*hbar^2*xl2*xl1^3)*x(i)+(-12*hbar^2*...
    xl1^5+12*hbar^2*xl2*xl1^4)))*xr2+(((720*m*xl1-720*m*xl2)*x(i)^4+(-2880*m*...
    xl1^2+2880*m*xl2*xl1)*x(i)^3+(4320*m*xl1^3-4320*m*xl2*xl1^2)*x(i)^2+(-2880*m*...
    xl1^4+2880*m*xl2*xl1^3-4*hbar^2)*x(i)+(720*m*xl1^5-720*m*xl2*xl1^4+16*hbar^2*...
    xl1-12*hbar^2*xl2))*xr1^5+((-2880*m*xl1+2880*m*xl2)*x(i)^5+(11520*m*xl1^2-11520*...
    m*xl2*xl1)*x(i)^4+(-17280*m*xl1^3+17280*m*xl2*xl1^2)*x(i)^3+(11520*m*...
    xl1^4-11520*m*xl2*xl1^3+16*hbar^2)*x(i)^2+(-2880*m*xl1^5+2880*m*xl2*xl1^4-68*...
    hbar^2*xl1+52*hbar^2*xl2)*x(i)+(4*hbar^2*xl1^2-4*hbar^2*xl2*xl1))*xr1^4+((4320*...
    m*xl1-4320*m*xl2)*x(i)^6+(-17280*m*xl1^2+17280*m*xl2*xl1)*x(i)^5+(25920*m*...
    xl1^3-25920*m*xl2*xl1^2)*x(i)^4+(-17280*m*xl1^4+17280*m*xl2*xl1^3-24*hbar^2)*...
    x(i)^3+(4320*m*xl1^5-4320*m*xl2*xl1^4+112*hbar^2*xl1-88*hbar^2*xl2)*x(i)^2+(-20*...
    hbar^2*xl1^2+20*hbar^2*xl2*xl1)*x(i)+(4*hbar^2*xl1^3-4*hbar^2*xl2*xl1^2))*...
    xr1^3+((-2880*m*xl1+2880*m*xl2)*x(i)^7+(11520*m*xl1^2-11520*m*xl2*xl1)*...
    x(i)^6+(-17280*m*xl1^3+17280*m*xl2*xl1^2)*x(i)^5+(11520*m*xl1^4-11520*m*xl2*...
    xl1^3+16*hbar^2)*x(i)^4+(-2880*m*xl1^5+2880*m*xl2*xl1^4-88*hbar^2*xl1+72*hbar^2*...
    xl2)*x(i)^3+(40*hbar^2*xl1^2-40*hbar^2*xl2*xl1)*x(i)^2+(-20*hbar^2*xl1^3+20*...
    hbar^2*xl2*xl1^2)*x(i)+(4*hbar^2*xl1^4-4*hbar^2*xl2*xl1^3))*xr1^2+((720*m*...
    xl1-720*m*xl2)*x(i)^8+(-2880*m*xl1^2+2880*m*xl2*xl1)*x(i)^7+(4320*m*xl1^3-4320*...
    m*xl2*xl1^2)*x(i)^6+(-2880*m*xl1^4+2880*m*xl2*xl1^3-4*hbar^2)*x(i)^5+(720*m*...
    xl1^5-720*m*xl2*xl1^4+44*hbar^2*xl1-40*hbar^2*xl2)*x(i)^4+(-88*hbar^2*xl1^2+88*...
    hbar^2*xl2*xl1)*x(i)^3+(112*hbar^2*xl1^3-112*hbar^2*xl2*xl1^2)*x(i)^2+(-68*...
    hbar^2*xl1^4+68*hbar^2*xl2*xl1^3)*x(i)+(16*hbar^2*xl1^5-16*hbar^2*xl2*xl1^4))*...
    xr1+((-4*hbar^2*xl1+4*hbar^2*xl2)*x(i)^5+(16*hbar^2*xl1^2-16*hbar^2*xl2*xl1)*...
    x(i)^4+(-24*hbar^2*xl1^3+24*hbar^2*xl2*xl1^2)*x(i)^3+(16*hbar^2*xl1^4-16*hbar^2*...
    xl2*xl1^3)*x(i)^2+(-4*hbar^2*xl1^5+4*hbar^2*xl2*xl1^4)*x(i))))*exp(5*x(i)-5*...
    c)^2+(((hbar^2*x(i)+(-4*hbar^2*xl1+3*hbar^2*xl2))*xr1^4+(-4*hbar^2*x(i)^2+(17*...
    hbar^2*xl1-13*hbar^2*xl2)*x(i)+(-1*hbar^2*xl1^2+hbar^2*xl2*xl1))*xr1^3+(6*...
    hbar^2*x(i)^3+(-28*hbar^2*xl1+22*hbar^2*xl2)*x(i)^2+(5*hbar^2*xl1^2-5*hbar^2*...
    xl2*xl1)*x(i)+(-1*hbar^2*xl1^3+hbar^2*xl2*xl1^2))*xr1^2+(-4*hbar^2*x(i)^4+(22*...
    hbar^2*xl1-18*hbar^2*xl2)*x(i)^3+(-10*hbar^2*xl1^2+10*hbar^2*xl2*xl1)*x(i)^2+(5*...
    hbar^2*xl1^3-5*hbar^2*xl2*xl1^2)*x(i)+(-1*hbar^2*xl1^4+hbar^2*xl2*xl1^3))*...
    xr1+(hbar^2*x(i)^5+(-10*hbar^2*xl1+9*hbar^2*xl2)*x(i)^4+(18*hbar^2*xl1^2-18*...
    hbar^2*xl2*xl1)*x(i)^3+(-22*hbar^2*xl1^3+22*hbar^2*xl2*xl1^2)*x(i)^2+(13*hbar^2*...
    xl1^4-13*hbar^2*xl2*xl1^3)*x(i)+(-3*hbar^2*xl1^5+3*hbar^2*xl2*xl1^4)))*xr2+((-1*...
    hbar^2*x(i)+(4*hbar^2*xl1-3*hbar^2*xl2))*xr1^5+(4*hbar^2*x(i)^2+(-17*hbar^2*...
    xl1+13*hbar^2*xl2)*x(i)+(hbar^2*xl1^2-1*hbar^2*xl2*xl1))*xr1^4+(-6*hbar^2*...
    x(i)^3+(28*hbar^2*xl1-22*hbar^2*xl2)*x(i)^2+(-5*hbar^2*xl1^2+5*hbar^2*xl2*xl1)*...
    x(i)+(hbar^2*xl1^3-1*hbar^2*xl2*xl1^2))*xr1^3+(4*hbar^2*x(i)^4+(-22*hbar^2*...
    xl1+18*hbar^2*xl2)*x(i)^3+(10*hbar^2*xl1^2-10*hbar^2*xl2*xl1)*x(i)^2+(-5*hbar^2*...
    xl1^3+5*hbar^2*xl2*xl1^2)*x(i)+(hbar^2*xl1^4-1*hbar^2*xl2*xl1^3))*xr1^2+(-1*...
    hbar^2*x(i)^5+(11*hbar^2*xl1-10*hbar^2*xl2)*x(i)^4+(-22*hbar^2*xl1^2+22*hbar^2*...
    xl2*xl1)*x(i)^3+(28*hbar^2*xl1^3-28*hbar^2*xl2*xl1^2)*x(i)^2+(-17*hbar^2*...
    xl1^4+17*hbar^2*xl2*xl1^3)*x(i)+(4*hbar^2*xl1^5-4*hbar^2*xl2*xl1^4))*xr1+((-1*...
    hbar^2*xl1+hbar^2*xl2)*x(i)^5+(4*hbar^2*xl1^2-4*hbar^2*xl2*xl1)*x(i)^4+(-6*...
    hbar^2*xl1^3+6*hbar^2*xl2*xl1^2)*x(i)^3+(4*hbar^2*xl1^4-4*hbar^2*xl2*xl1^3)*...
    x(i)^2+(-1*hbar^2*xl1^5+hbar^2*xl2*xl1^4)*x(i)))))/(((((2*m^2*xl1-2*m^2*xl2)*...
    x(i)^4+(-8*m^2*xl1^2+8*m^2*xl2*xl1)*x(i)^3+(12*m^2*xl1^3-12*m^2*xl2*xl1^2)*...
    x(i)^2+(-8*m^2*xl1^4+8*m^2*xl2*xl1^3)*x(i)+(2*m^2*xl1^5-2*m^2*xl2*xl1^4))*...
    xr1^4+((-8*m^2*xl1+8*m^2*xl2)*x(i)^5+(32*m^2*xl1^2-32*m^2*xl2*xl1)*x(i)^4+(-48*...
    m^2*xl1^3+48*m^2*xl2*xl1^2)*x(i)^3+(32*m^2*xl1^4-32*m^2*xl2*xl1^3)*x(i)^2+(-8*...
    m^2*xl1^5+8*m^2*xl2*xl1^4)*x(i))*xr1^3+((12*m^2*xl1-12*m^2*xl2)*x(i)^6+(-48*m^2*...
    xl1^2+48*m^2*xl2*xl1)*x(i)^5+(72*m^2*xl1^3-72*m^2*xl2*xl1^2)*x(i)^4+(-48*m^2*...
    xl1^4+48*m^2*xl2*xl1^3)*x(i)^3+(12*m^2*xl1^5-12*m^2*xl2*xl1^4)*x(i)^2)*...
    xr1^2+((-8*m^2*xl1+8*m^2*xl2)*x(i)^7+(32*m^2*xl1^2-32*m^2*xl2*xl1)*x(i)^6+(-48*...
    m^2*xl1^3+48*m^2*xl2*xl1^2)*x(i)^5+(32*m^2*xl1^4-32*m^2*xl2*xl1^3)*x(i)^4+(-8*...
    m^2*xl1^5+8*m^2*xl2*xl1^4)*x(i)^3)*xr1+((2*m^2*xl1-2*m^2*xl2)*x(i)^8+(-8*m^2*...
    xl1^2+8*m^2*xl2*xl1)*x(i)^7+(12*m^2*xl1^3-12*m^2*xl2*xl1^2)*x(i)^6+(-8*m^2*...
    xl1^4+8*m^2*xl2*xl1^3)*x(i)^5+(2*m^2*xl1^5-2*m^2*xl2*xl1^4)*x(i)^4))*xr2+(((-2*...
    m^2*xl1+2*m^2*xl2)*x(i)^4+(8*m^2*xl1^2-8*m^2*xl2*xl1)*x(i)^3+(-12*m^2*xl1^3+12*...
    m^2*xl2*xl1^2)*x(i)^2+(8*m^2*xl1^4-8*m^2*xl2*xl1^3)*x(i)+(-2*m^2*xl1^5+2*m^2*...
    xl2*xl1^4))*xr1^5+((8*m^2*xl1-8*m^2*xl2)*x(i)^5+(-32*m^2*xl1^2+32*m^2*xl2*xl1)*...
    x(i)^4+(48*m^2*xl1^3-48*m^2*xl2*xl1^2)*x(i)^3+(-32*m^2*xl1^4+32*m^2*xl2*xl1^3)*...
    x(i)^2+(8*m^2*xl1^5-8*m^2*xl2*xl1^4)*x(i))*xr1^4+((-12*m^2*xl1+12*m^2*xl2)*...
    x(i)^6+(48*m^2*xl1^2-48*m^2*xl2*xl1)*x(i)^5+(-72*m^2*xl1^3+72*m^2*xl2*xl1^2)*...
    x(i)^4+(48*m^2*xl1^4-48*m^2*xl2*xl1^3)*x(i)^3+(-12*m^2*xl1^5+12*m^2*xl2*xl1^4)*...
    x(i)^2)*xr1^3+((8*m^2*xl1-8*m^2*xl2)*x(i)^7+(-32*m^2*xl1^2+32*m^2*xl2*xl1)*...
    x(i)^6+(48*m^2*xl1^3-48*m^2*xl2*xl1^2)*x(i)^5+(-32*m^2*xl1^4+32*m^2*xl2*xl1^3)*...
    x(i)^4+(8*m^2*xl1^5-8*m^2*xl2*xl1^4)*x(i)^3)*xr1^2+((-2*m^2*xl1+2*m^2*xl2)*...
    x(i)^8+(8*m^2*xl1^2-8*m^2*xl2*xl1)*x(i)^7+(-12*m^2*xl1^3+12*m^2*xl2*xl1^2)*...
    x(i)^6+(8*m^2*xl1^4-8*m^2*xl2*xl1^3)*x(i)^5+(-2*m^2*xl1^5+2*m^2*xl2*xl1^4)*...
    x(i)^4)*xr1))*exp(5*x(i)-5*c)^8+((((8*m^2*xl1-8*m^2*xl2)*x(i)^4+(-32*m^2*...
    xl1^2+32*m^2*xl2*xl1)*x(i)^3+(48*m^2*xl1^3-48*m^2*xl2*xl1^2)*x(i)^2+(-32*m^2*...
    xl1^4+32*m^2*xl2*xl1^3)*x(i)+(8*m^2*xl1^5-8*m^2*xl2*xl1^4))*xr1^4+((-32*m^2*...
    xl1+32*m^2*xl2)*x(i)^5+(128*m^2*xl1^2-128*m^2*xl2*xl1)*x(i)^4+(-192*m^2*...
    xl1^3+192*m^2*xl2*xl1^2)*x(i)^3+(128*m^2*xl1^4-128*m^2*xl2*xl1^3)*x(i)^2+(-32*...
    m^2*xl1^5+32*m^2*xl2*xl1^4)*x(i))*xr1^3+((48*m^2*xl1-48*m^2*xl2)*x(i)^6+(-192*...
    m^2*xl1^2+192*m^2*xl2*xl1)*x(i)^5+(288*m^2*xl1^3-288*m^2*xl2*xl1^2)*...
    x(i)^4+(-192*m^2*xl1^4+192*m^2*xl2*xl1^3)*x(i)^3+(48*m^2*xl1^5-48*m^2*xl2*...
    xl1^4)*x(i)^2)*xr1^2+((-32*m^2*xl1+32*m^2*xl2)*x(i)^7+(128*m^2*xl1^2-128*m^2*...
    xl2*xl1)*x(i)^6+(-192*m^2*xl1^3+192*m^2*xl2*xl1^2)*x(i)^5+(128*m^2*xl1^4-128*...
    m^2*xl2*xl1^3)*x(i)^4+(-32*m^2*xl1^5+32*m^2*xl2*xl1^4)*x(i)^3)*xr1+((8*m^2*...
    xl1-8*m^2*xl2)*x(i)^8+(-32*m^2*xl1^2+32*m^2*xl2*xl1)*x(i)^7+(48*m^2*xl1^3-48*...
    m^2*xl2*xl1^2)*x(i)^6+(-32*m^2*xl1^4+32*m^2*xl2*xl1^3)*x(i)^5+(8*m^2*xl1^5-8*...
    m^2*xl2*xl1^4)*x(i)^4))*xr2+(((-8*m^2*xl1+8*m^2*xl2)*x(i)^4+(32*m^2*xl1^2-32*...
    m^2*xl2*xl1)*x(i)^3+(-48*m^2*xl1^3+48*m^2*xl2*xl1^2)*x(i)^2+(32*m^2*xl1^4-32*...
    m^2*xl2*xl1^3)*x(i)+(-8*m^2*xl1^5+8*m^2*xl2*xl1^4))*xr1^5+((32*m^2*xl1-32*m^2*...
    xl2)*x(i)^5+(-128*m^2*xl1^2+128*m^2*xl2*xl1)*x(i)^4+(192*m^2*xl1^3-192*m^2*xl2*...
    xl1^2)*x(i)^3+(-128*m^2*xl1^4+128*m^2*xl2*xl1^3)*x(i)^2+(32*m^2*xl1^5-32*m^2*...
    xl2*xl1^4)*x(i))*xr1^4+((-48*m^2*xl1+48*m^2*xl2)*x(i)^6+(192*m^2*xl1^2-192*m^2*...
    xl2*xl1)*x(i)^5+(-288*m^2*xl1^3+288*m^2*xl2*xl1^2)*x(i)^4+(192*m^2*xl1^4-192*...
    m^2*xl2*xl1^3)*x(i)^3+(-48*m^2*xl1^5+48*m^2*xl2*xl1^4)*x(i)^2)*xr1^3+((32*m^2*...
    xl1-32*m^2*xl2)*x(i)^7+(-128*m^2*xl1^2+128*m^2*xl2*xl1)*x(i)^6+(192*m^2*...
    xl1^3-192*m^2*xl2*xl1^2)*x(i)^5+(-128*m^2*xl1^4+128*m^2*xl2*xl1^3)*x(i)^4+(32*...
    m^2*xl1^5-32*m^2*xl2*xl1^4)*x(i)^3)*xr1^2+((-8*m^2*xl1+8*m^2*xl2)*x(i)^8+(32*...
    m^2*xl1^2-32*m^2*xl2*xl1)*x(i)^7+(-48*m^2*xl1^3+48*m^2*xl2*xl1^2)*x(i)^6+(32*...
    m^2*xl1^4-32*m^2*xl2*xl1^3)*x(i)^5+(-8*m^2*xl1^5+8*m^2*xl2*xl1^4)*x(i)^4)*xr1))*...
    exp(5*x(i)-5*c)^6+((((12*m^2*xl1-12*m^2*xl2)*x(i)^4+(-48*m^2*xl1^2+48*m^2*xl2*...
    xl1)*x(i)^3+(72*m^2*xl1^3-72*m^2*xl2*xl1^2)*x(i)^2+(-48*m^2*xl1^4+48*m^2*xl2*...
    xl1^3)*x(i)+(12*m^2*xl1^5-12*m^2*xl2*xl1^4))*xr1^4+((-48*m^2*xl1+48*m^2*xl2)*...
    x(i)^5+(192*m^2*xl1^2-192*m^2*xl2*xl1)*x(i)^4+(-288*m^2*xl1^3+288*m^2*xl2*...
    xl1^2)*x(i)^3+(192*m^2*xl1^4-192*m^2*xl2*xl1^3)*x(i)^2+(-48*m^2*xl1^5+48*m^2*...
    xl2*xl1^4)*x(i))*xr1^3+((72*m^2*xl1-72*m^2*xl2)*x(i)^6+(-288*m^2*xl1^2+288*m^2*...
    xl2*xl1)*x(i)^5+(432*m^2*xl1^3-432*m^2*xl2*xl1^2)*x(i)^4+(-288*m^2*xl1^4+288*...
    m^2*xl2*xl1^3)*x(i)^3+(72*m^2*xl1^5-72*m^2*xl2*xl1^4)*x(i)^2)*xr1^2+((-48*m^2*...
    xl1+48*m^2*xl2)*x(i)^7+(192*m^2*xl1^2-192*m^2*xl2*xl1)*x(i)^6+(-288*m^2*...
    xl1^3+288*m^2*xl2*xl1^2)*x(i)^5+(192*m^2*xl1^4-192*m^2*xl2*xl1^3)*x(i)^4+(-48*...
    m^2*xl1^5+48*m^2*xl2*xl1^4)*x(i)^3)*xr1+((12*m^2*xl1-12*m^2*xl2)*x(i)^8+(-48*...
    m^2*xl1^2+48*m^2*xl2*xl1)*x(i)^7+(72*m^2*xl1^3-72*m^2*xl2*xl1^2)*x(i)^6+(-48*...
    m^2*xl1^4+48*m^2*xl2*xl1^3)*x(i)^5+(12*m^2*xl1^5-12*m^2*xl2*xl1^4)*x(i)^4))*...
    xr2+(((-12*m^2*xl1+12*m^2*xl2)*x(i)^4+(48*m^2*xl1^2-48*m^2*xl2*xl1)*x(i)^3+(-72*...
    m^2*xl1^3+72*m^2*xl2*xl1^2)*x(i)^2+(48*m^2*xl1^4-48*m^2*xl2*xl1^3)*x(i)+(-12*...
    m^2*xl1^5+12*m^2*xl2*xl1^4))*xr1^5+((48*m^2*xl1-48*m^2*xl2)*x(i)^5+(-192*m^2*...
    xl1^2+192*m^2*xl2*xl1)*x(i)^4+(288*m^2*xl1^3-288*m^2*xl2*xl1^2)*x(i)^3+(-192*...
    m^2*xl1^4+192*m^2*xl2*xl1^3)*x(i)^2+(48*m^2*xl1^5-48*m^2*xl2*xl1^4)*x(i))*...
    xr1^4+((-72*m^2*xl1+72*m^2*xl2)*x(i)^6+(288*m^2*xl1^2-288*m^2*xl2*xl1)*...
    x(i)^5+(-432*m^2*xl1^3+432*m^2*xl2*xl1^2)*x(i)^4+(288*m^2*xl1^4-288*m^2*xl2*...
    xl1^3)*x(i)^3+(-72*m^2*xl1^5+72*m^2*xl2*xl1^4)*x(i)^2)*xr1^3+((48*m^2*xl1-48*...
    m^2*xl2)*x(i)^7+(-192*m^2*xl1^2+192*m^2*xl2*xl1)*x(i)^6+(288*m^2*xl1^3-288*m^2*...
    xl2*xl1^2)*x(i)^5+(-192*m^2*xl1^4+192*m^2*xl2*xl1^3)*x(i)^4+(48*m^2*xl1^5-48*...
    m^2*xl2*xl1^4)*x(i)^3)*xr1^2+((-12*m^2*xl1+12*m^2*xl2)*x(i)^8+(48*m^2*xl1^2-48*...
    m^2*xl2*xl1)*x(i)^7+(-72*m^2*xl1^3+72*m^2*xl2*xl1^2)*x(i)^6+(48*m^2*xl1^4-48*...
    m^2*xl2*xl1^3)*x(i)^5+(-12*m^2*xl1^5+12*m^2*xl2*xl1^4)*x(i)^4)*xr1))*exp(5*...
    x(i)-5*c)^4+((((8*m^2*xl1-8*m^2*xl2)*x(i)^4+(-32*m^2*xl1^2+32*m^2*xl2*xl1)*...
    x(i)^3+(48*m^2*xl1^3-48*m^2*xl2*xl1^2)*x(i)^2+(-32*m^2*xl1^4+32*m^2*xl2*xl1^3)*...
    x(i)+(8*m^2*xl1^5-8*m^2*xl2*xl1^4))*xr1^4+((-32*m^2*xl1+32*m^2*xl2)*x(i)^5+(128*...
    m^2*xl1^2-128*m^2*xl2*xl1)*x(i)^4+(-192*m^2*xl1^3+192*m^2*xl2*xl1^2)*...
    x(i)^3+(128*m^2*xl1^4-128*m^2*xl2*xl1^3)*x(i)^2+(-32*m^2*xl1^5+32*m^2*xl2*...
    xl1^4)*x(i))*xr1^3+((48*m^2*xl1-48*m^2*xl2)*x(i)^6+(-192*m^2*xl1^2+192*m^2*xl2*...
    xl1)*x(i)^5+(288*m^2*xl1^3-288*m^2*xl2*xl1^2)*x(i)^4+(-192*m^2*xl1^4+192*m^2*...
    xl2*xl1^3)*x(i)^3+(48*m^2*xl1^5-48*m^2*xl2*xl1^4)*x(i)^2)*xr1^2+((-32*m^2*...
    xl1+32*m^2*xl2)*x(i)^7+(128*m^2*xl1^2-128*m^2*xl2*xl1)*x(i)^6+(-192*m^2*...
    xl1^3+192*m^2*xl2*xl1^2)*x(i)^5+(128*m^2*xl1^4-128*m^2*xl2*xl1^3)*x(i)^4+(-32*...
    m^2*xl1^5+32*m^2*xl2*xl1^4)*x(i)^3)*xr1+((8*m^2*xl1-8*m^2*xl2)*x(i)^8+(-32*m^2*...
    xl1^2+32*m^2*xl2*xl1)*x(i)^7+(48*m^2*xl1^3-48*m^2*xl2*xl1^2)*x(i)^6+(-32*m^2*...
    xl1^4+32*m^2*xl2*xl1^3)*x(i)^5+(8*m^2*xl1^5-8*m^2*xl2*xl1^4)*x(i)^4))*xr2+(((-8*...
    m^2*xl1+8*m^2*xl2)*x(i)^4+(32*m^2*xl1^2-32*m^2*xl2*xl1)*x(i)^3+(-48*m^2*...
    xl1^3+48*m^2*xl2*xl1^2)*x(i)^2+(32*m^2*xl1^4-32*m^2*xl2*xl1^3)*x(i)+(-8*m^2*...
    xl1^5+8*m^2*xl2*xl1^4))*xr1^5+((32*m^2*xl1-32*m^2*xl2)*x(i)^5+(-128*m^2*...
    xl1^2+128*m^2*xl2*xl1)*x(i)^4+(192*m^2*xl1^3-192*m^2*xl2*xl1^2)*x(i)^3+(-128*...
    m^2*xl1^4+128*m^2*xl2*xl1^3)*x(i)^2+(32*m^2*xl1^5-32*m^2*xl2*xl1^4)*x(i))*...
    xr1^4+((-48*m^2*xl1+48*m^2*xl2)*x(i)^6+(192*m^2*xl1^2-192*m^2*xl2*xl1)*...
    x(i)^5+(-288*m^2*xl1^3+288*m^2*xl2*xl1^2)*x(i)^4+(192*m^2*xl1^4-192*m^2*xl2*...
    xl1^3)*x(i)^3+(-48*m^2*xl1^5+48*m^2*xl2*xl1^4)*x(i)^2)*xr1^3+((32*m^2*xl1-32*...
    m^2*xl2)*x(i)^7+(-128*m^2*xl1^2+128*m^2*xl2*xl1)*x(i)^6+(192*m^2*xl1^3-192*m^2*...
    xl2*xl1^2)*x(i)^5+(-128*m^2*xl1^4+128*m^2*xl2*xl1^3)*x(i)^4+(32*m^2*xl1^5-32*...
    m^2*xl2*xl1^4)*x(i)^3)*xr1^2+((-8*m^2*xl1+8*m^2*xl2)*x(i)^8+(32*m^2*xl1^2-32*...
    m^2*xl2*xl1)*x(i)^7+(-48*m^2*xl1^3+48*m^2*xl2*xl1^2)*x(i)^6+(32*m^2*xl1^4-32*...
    m^2*xl2*xl1^3)*x(i)^5+(-8*m^2*xl1^5+8*m^2*xl2*xl1^4)*x(i)^4)*xr1))*exp(5*x(i)-5*...
    c)^2+((((2*m^2*xl1-2*m^2*xl2)*x(i)^4+(-8*m^2*xl1^2+8*m^2*xl2*xl1)*x(i)^3+(12*...
    m^2*xl1^3-12*m^2*xl2*xl1^2)*x(i)^2+(-8*m^2*xl1^4+8*m^2*xl2*xl1^3)*x(i)+(2*m^2*...
    xl1^5-2*m^2*xl2*xl1^4))*xr1^4+((-8*m^2*xl1+8*m^2*xl2)*x(i)^5+(32*m^2*xl1^2-32*...
    m^2*xl2*xl1)*x(i)^4+(-48*m^2*xl1^3+48*m^2*xl2*xl1^2)*x(i)^3+(32*m^2*xl1^4-32*...
    m^2*xl2*xl1^3)*x(i)^2+(-8*m^2*xl1^5+8*m^2*xl2*xl1^4)*x(i))*xr1^3+((12*m^2*...
    xl1-12*m^2*xl2)*x(i)^6+(-48*m^2*xl1^2+48*m^2*xl2*xl1)*x(i)^5+(72*m^2*xl1^3-72*...
    m^2*xl2*xl1^2)*x(i)^4+(-48*m^2*xl1^4+48*m^2*xl2*xl1^3)*x(i)^3+(12*m^2*xl1^5-12*...
    m^2*xl2*xl1^4)*x(i)^2)*xr1^2+((-8*m^2*xl1+8*m^2*xl2)*x(i)^7+(32*m^2*xl1^2-32*...
    m^2*xl2*xl1)*x(i)^6+(-48*m^2*xl1^3+48*m^2*xl2*xl1^2)*x(i)^5+(32*m^2*xl1^4-32*...
    m^2*xl2*xl1^3)*x(i)^4+(-8*m^2*xl1^5+8*m^2*xl2*xl1^4)*x(i)^3)*xr1+((2*m^2*xl1-2*...
    m^2*xl2)*x(i)^8+(-8*m^2*xl1^2+8*m^2*xl2*xl1)*x(i)^7+(12*m^2*xl1^3-12*m^2*xl2*...
    xl1^2)*x(i)^6+(-8*m^2*xl1^4+8*m^2*xl2*xl1^3)*x(i)^5+(2*m^2*xl1^5-2*m^2*xl2*...
    xl1^4)*x(i)^4))*xr2+(((-2*m^2*xl1+2*m^2*xl2)*x(i)^4+(8*m^2*xl1^2-8*m^2*xl2*xl1)*...
    x(i)^3+(-12*m^2*xl1^3+12*m^2*xl2*xl1^2)*x(i)^2+(8*m^2*xl1^4-8*m^2*xl2*xl1^3)*...
    x(i)+(-2*m^2*xl1^5+2*m^2*xl2*xl1^4))*xr1^5+((8*m^2*xl1-8*m^2*xl2)*x(i)^5+(-32*...
    m^2*xl1^2+32*m^2*xl2*xl1)*x(i)^4+(48*m^2*xl1^3-48*m^2*xl2*xl1^2)*x(i)^3+(-32*...
    m^2*xl1^4+32*m^2*xl2*xl1^3)*x(i)^2+(8*m^2*xl1^5-8*m^2*xl2*xl1^4)*x(i))*...
    xr1^4+((-12*m^2*xl1+12*m^2*xl2)*x(i)^6+(48*m^2*xl1^2-48*m^2*xl2*xl1)*...
    x(i)^5+(-72*m^2*xl1^3+72*m^2*xl2*xl1^2)*x(i)^4+(48*m^2*xl1^4-48*m^2*xl2*xl1^3)*...
    x(i)^3+(-12*m^2*xl1^5+12*m^2*xl2*xl1^4)*x(i)^2)*xr1^3+((8*m^2*xl1-8*m^2*xl2)*...
    x(i)^7+(-32*m^2*xl1^2+32*m^2*xl2*xl1)*x(i)^6+(48*m^2*xl1^3-48*m^2*xl2*xl1^2)*...
    x(i)^5+(-32*m^2*xl1^4+32*m^2*xl2*xl1^3)*x(i)^4+(8*m^2*xl1^5-8*m^2*xl2*xl1^4)*...
    x(i)^3)*xr1^2+((-2*m^2*xl1+2*m^2*xl2)*x(i)^8+(8*m^2*xl1^2-8*m^2*xl2*xl1)*...
    x(i)^7+(-12*m^2*xl1^3+12*m^2*xl2*xl1^2)*x(i)^6+(8*m^2*xl1^4-8*m^2*xl2*xl1^3)*...
    x(i)^5+(-2*m^2*xl1^5+2*m^2*xl2*xl1^4)*x(i)^4)*xr1)));
endfunction

%sage print fricas2octave('Jp1','jac1(i,1)')

function result=Jp1(x,i)
  global l2 l1 p1 p2;
  global N m hbar a b c;
  if (i-2<1)
    xl2=l2;
  else
    xl2=x(i-2);
  endif;
  if (i-1<1)
    xl1=l1;
  else
    xl1=x(i-1);
  endif;
  if (i+2>N)
    xr2=p2;
  else
    xr2=x(i+2);
  endif;
  if (i+1>N)
    xr1=p1;
  else
    xr1=x(i+1);
  endif;
  result=((hbar^2*xr1^2+(-4*hbar^2*x(i)+2*hbar^2*xl1)*xr1+(9*hbar^2*x(i)^2-14*hbar^2*xl1*...
    x(i)+6*hbar^2*xl1^2))*xr2^2+(-2*hbar^2*xr1^3+(8*hbar^2*x(i)-4*hbar^2*xl1)*...
    xr1^2+(-20*hbar^2*x(i)^2+32*hbar^2*xl1*x(i)-14*hbar^2*xl1^2)*xr1+(2*hbar^2*...
    x(i)^3-4*hbar^2*xl1*x(i)^2+2*hbar^2*xl1^2*x(i)))*xr2+(hbar^2*xr1^4+(-4*hbar^2*...
    x(i)+2*hbar^2*xl1)*xr1^3+(12*hbar^2*x(i)^2-20*hbar^2*xl1*x(i)+9*hbar^2*xl1^2)*...
    xr1^2+(-4*hbar^2*x(i)^3+8*hbar^2*xl1*x(i)^2-4*hbar^2*xl1^2*x(i))*xr1+(hbar^2*...
    x(i)^4-2*hbar^2*xl1*x(i)^3+hbar^2*xl1^2*x(i)^2)))/(((4*m^2*x(i)^2-8*m^2*xl1*...
    x(i)+4*m^2*xl1^2)*xr1^4+(-16*m^2*x(i)^3+32*m^2*xl1*x(i)^2-16*m^2*xl1^2*x(i))*...
    xr1^3+(24*m^2*x(i)^4-48*m^2*xl1*x(i)^3+24*m^2*xl1^2*x(i)^2)*xr1^2+(-16*m^2*...
    x(i)^5+32*m^2*xl1*x(i)^4-16*m^2*xl1^2*x(i)^3)*xr1+(4*m^2*x(i)^6-8*m^2*xl1*...
    x(i)^5+4*m^2*xl1^2*x(i)^4))*xr2^2+((-8*m^2*x(i)^2+16*m^2*xl1*x(i)-8*m^2*xl1^2)*...
    xr1^5+(32*m^2*x(i)^3-64*m^2*xl1*x(i)^2+32*m^2*xl1^2*x(i))*xr1^4+(-48*m^2*...
    x(i)^4+96*m^2*xl1*x(i)^3-48*m^2*xl1^2*x(i)^2)*xr1^3+(32*m^2*x(i)^5-64*m^2*xl1*...
    x(i)^4+32*m^2*xl1^2*x(i)^3)*xr1^2+(-8*m^2*x(i)^6+16*m^2*xl1*x(i)^5-8*m^2*xl1^2*...
    x(i)^4)*xr1)*xr2+((4*m^2*x(i)^2-8*m^2*xl1*x(i)+4*m^2*xl1^2)*xr1^6+(-16*m^2*...
    x(i)^3+32*m^2*xl1*x(i)^2-16*m^2*xl1^2*x(i))*xr1^5+(24*m^2*x(i)^4-48*m^2*xl1*...
    x(i)^3+24*m^2*xl1^2*x(i)^2)*xr1^4+(-16*m^2*x(i)^5+32*m^2*xl1*x(i)^4-16*m^2*...
    xl1^2*x(i)^3)*xr1^3+(4*m^2*x(i)^6-8*m^2*xl1*x(i)^5+4*m^2*xl1^2*x(i)^4)*xr1^2));
endfunction

%sage print fricas2octave('Jp2','jac1(i,2)')

function result=Jp2(x,i)
  global l2 l1 p1 p2;
  global N m hbar a b c;
  if (i-2<1)
    xl2=l2;
  else
    xl2=x(i-2);
  endif;
  if (i-1<1)
    xl1=l1;
  else
    xl1=x(i-1);
  endif;
  if (i+2>N)
    xr2=p2;
  else
    xr2=x(i+2);
  endif;
  if (i+1>N)
    xr1=p1;
  else
    xr1=x(i+1);
  endif;
  result=(-1*hbar^2)/((4*m^2*xr1^2-8*m^2*x(i)*xr1+4*m^2*x(i)^2)*xr2^2+(-8*m^2*xr1^3+16*...
    m^2*x(i)*xr1^2-8*m^2*x(i)^2*xr1)*xr2+(4*m^2*xr1^4-8*m^2*x(i)*xr1^3+4*m^2*x(i)^2*...
    xr1^2));
endfunction

function result=jac1(x,i,j)
  global m hbar a;
  if (j==-2)
    result=Jm2(x,i);
  elseif (j==-1)
    result=Jm1(x,i);
  elseif (j==0)
    result=J0(x,i);
  elseif (j==1)
    result=Jp1(x,i);
  elseif (j==2)
    result=Jp2(x,i);
  else
    result=0;
  endif
endfunction

jac1([X0;Y0],3,2)

ans = -2.50144

function result = jac0(x,i,j)
  global N m hbar a;
  if (i==j-N)
    result = 1;
  elseif (i<=N)
    result = 0;
  elseif (i<N+j-2)
    result = 0;
  else
    result = jac1(x,i-N,i-N-j);
  endif
endfunction

jac0([X0;Y0],15,5)

ans = 0

function jac=jacm(x,t)
  global N
  jac=zeros(2*N,2*N);
  for i=1:2*N
    for j=1:2*N
      jac(i,j)=jac0(x,i,j);
    endfor
  endfor
endfunction

Perform the time updates

t0 = clock ();
trec = linspace (0, 10, 240)';

#Xrec = lsode (@deq, [X0;Y0], trec)';
Xrec = lsode ({@deq,@jacm}, [X0;Y0], trec)';
elapsed_time = etime (clock (), t0)
T=size(trec,1);
lsode_options()

elapsed_time = 108.378

Options for LSODE include:

  keyword                                             value
  -------                                             -----
  absolute tolerance                                  1.49012e-08
  relative tolerance                                  1.49012e-08
  integration method                                  stiff
  initial step size                                   -1
  maximum order                                       -1
  maximum step size                                   -1
  minimum step size                                   0
  step limit                                          100000

h = figure('visible', 'off');
plot(trec,Xrec(1:N,1:T));
saveas(h,"figure2.png");
%ylim([-2,2]);

%sage salvus.file("figure2.png")

h = figure('visible', 'off');
%Plot the density at the end of the simulation
Dist=((Xrec(2:N,T)-Xrec(1:(N-1),T)));
Prob=(1/N)*Dist.^(-1);
plot(Xrec(2:N,T),Prob);
saveas(h,"figure3.png");
%xlim([-1.5,1.5]); ylim([0,1]);

%sage salvus.file("figure3.png")

h = figure('visible', 'off');
%Plot the density at 3/4 through the simulation
Tsamp=3*T/4;
Dist=((Xrec(2:N,Tsamp)-Xrec(1:(N-1),Tsamp)));
Prob=(1/N)*Dist.^(-1);

plot(Xrec(2:N,Tsamp),Prob);
saveas(h,"figure4.png");
%xlim([-1.5,1.5]); ylim([0,1]);

%sage salvus.file("figure4.png")

h = figure('visible', 'off');
Dist=((Xrec(2:N,T)-Xrec(1:(N-1),T)));
Prob=(1/N)*Dist.^(-1);
plot(xvalues(2:N,1),Pin);
hold on
Dist=((Xrec(2:N,T)-Xrec(1:(N-1),T)));
Prob=(1/N)*Dist.^(-1);
plot(Xrec(2:N,T),Prob,"r");
saveas(h,"figure5.png");
hold off

%sage salvus.file("figure5.png")

figure("visible","off");
moviename="MIW1Dv4.mp4";
movie("init",moviename)
for Tsamp=1:10:T
Dist=((Xrec(2:N,Tsamp)-Xrec(1:(N-1),Tsamp)));
Prob=Dist.^(-1);
plot(Xrec(2:N,Tsamp),Prob);
xlim([-20,20]); ylim([0,20]);
movie("add",moviename);
endfor
movie("close",moviename,4); %last arg is frame rate
T

T = 240

%sage salvus.file("MIW1Dv4.mp4")